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Sommario 

La base di Shiff protonata 4-( cyclo-pent-2-enyldiene )-3,4-dihydro-2H-pyrrolinmm (CPP) è una 
molecola biomimetica proposta come modello per il fotorecettore del sistema visivo, e come base 
per lo sviluppo di una nuova classe di interruttori molecolari (switch) fotochimici. 

Nella presente tesi è stata studiata a livello teorico computazionale la fotoisomerizzazione del 
CPP con il metodo ab initio Multi-Configurational Self Consistent Field (MCSCF) nella modalità 
Complete Active Space (CASSCF). 

Questi calcoli hanno permesso di ottenere una mappatura delle Superfici di Energia Potenziale 
(PES) per stato eccitato (Si) e fondamentale (So) lungo quattro coordinate interne chimicamente 
significative. 

L'analisi di questi dati ha permesso la riduzione della dimensionalità del problema e limitare 
gli studi successivi a tre sole coordinate interne lungo le quali il sistema raggiunge l'intersezione 
conica tra i due stati elettronici (CI) responsabile del decadimento nonradiativo e della conseguente 
isomerizzazione del CPP. E emersa una correlazione tra la coordinata di torsione tra gli anelli, 
la piramidalizzazione dell'atomo di carbonio pirrolinico che li congiunge, ed il gradiente della 
PES per lo stato Sq nella regione contigua all'intersezione conica. Tale correlazione è confermata 
da un modello semplificato Extended Hùckel per il sistema tt e dalla morfologia degli orbitali 
appartenenti allo spazio attivo. 

La matrice del momento di dipolo elettrico relativa agli stati So ed Si è stata calcolata in tutti 
i punti delle superfici ed utilizzata per trasformare le PES adiabatiche, fortemente discontinue 
intorno alla intersezione, in PES diabatiche più continue rispetto alle coordinate nucleari. Per le 
due superfici diabatiche ed il loro accoppiamento è stata quindi ricercata una forma funzionale 
mediante fit lineare. 

Le superfici diabatiche ottenute sono state utilizzate per simulare la dinamica semiclassica della 
fotoisomerizzazione sulla scala dei picosecondi, investigando il meccanismo di reazione ed i fattori 
che influenzano la resa quantica. 
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1 Introduzione 



Questa Tesi riguarda il primo istante dell'interazione tra l'occhio umano e la luce: la scintilla 
fotochimica capace di saettare fin dentro al nostro cervello producendo la visione. Per indagare 
questo attimo così prezioso per la nostra vita abbiamo impiegato le più raffinate teorie di chimica 
quantistica, tecniche computazionali avanzate e cluster di computer per calcolo ad alte prestazioni. 

I nostri sforzi si sono concentrati su una parte piccolissima di ciò che si conosce del sistema visivo 
umano, cioè un modello biomimetico per il retinale, cromoforo responsabile dell'assorbimento 
luminoso: il CPP, 4-(cyclo-pent-2-enyldiene)-3,4-dihydro-2H-pyrrolinium. 

Si stima che il corpo umano contenga circa un miliardo di cellule per grammo, ciascuna composta 
da qualcosa come mille miliardi di atomi. I ventidue atomi di cui è composto il CPP possono allora 
sembrare insignificanti rispetto a sette grammi di occhio (7E+21 atomi). Inoltre con i metodi 
impiegati abbiamo potuto studiare una durata temporale di solo alcune centinaia di femtosecondi 
(10~^^5) da quando la luce colpisce il recettore: se paragonati ai cento millisecondi necessari alla 
sola attivazione della cellula fotosensibile, appaiono anch'essi trascurabili. 

L'utilità di uno studio altamente accurato su un modello così apparentemente limitato si rivela 
quando venga integrato con sforzi di ricerca più estesi e da esso sia possibile intuire principi più 
generali. Nella fattispecie i risultati conseguiti si inseriscono in un progetto di ricerca più ampio, 
volto allo studio dei processi fotochimici della proteina che costituisce il fotorecettore umano, la 
Rodopsina. 

Proprio in virtù della velocità e dell'efficienza del sistema biologico che modellizza, il CPP è 
stato proposto come scheletro per una nuova famiglia di interruttori e motori molecolari foto- 
chimici [42] ultraveloci, ovvero molecole capaci di compiere lavoro utile attraverso cambiamenti 
conformazionali quando stimolate da un impulso luminoso. Comprendere il comportamento del 
CPP dal punto di vista teorico è dunque cruciale per indirizzare la sintesi e l'ingegnerizzazione di 
nuovi motori. 

Le prossime sezioni introduttive procedono dall'occhio umano fino al livello molecolare per 
ingrandimenti successivi, ampliando poi il campo alle reazioni ultraveloci in genere (femtochimica) 
ed ai motori molecolari. Si innesta così il nostro lavoro in più di 100 anni di ricerca scientifica, 
permettendo al lettore di comprenderne gli obiettivi; lo stile sarà pertanto più narrativo che 
scientifico. Nel Capitolo 2 viene impostato il modello per il CPP, si descrive il protocollo di 
calcolo e vengono presentati i risultati preliminari. Il Capitolo 3 illustra le superfici di energia 
potenziale, la costruzione delle superfici diabatiche ed il fit lineare. Il capitolo 4 riporta i primi 
risultati di dinamica molecolare semiclassica per la fotoisomerizzazione. Nelle appendici sono 
descritti i programmi sviluppati per risolvere le difficoltà computazionali che abbiamo incontrato. 



I Introduzione 

1.1 II sistema visivo 

II complesso meccanismo della visione comincia nel momento in cui la luce colpisce la cornea 
del nostro occhio (fig.1.1), dove viene deviata dal diverso indice di rifrazione e focalizzata dalla 
curvatura della pellicola corneale. La luce attraversa poi l'umor acqueo, che riempie la cornea 
mantenendone la forma emisferica, e la pupilla che ne regola l'intensitàà. 




Figura 1.1: L'occhio umano 

Giungendo al cristallino, una sfera appiattita costituita da fasci sovrapposti di fibre trasparenti 
la cui forma è modulata dai muscoli ciliari. Il cristallino funziona da lente biconvessa con curvature 
differenti e converge i raggi sulla retina, ove si compone una immagine rovesciata di ciò che l'occhio 
sta osservando. L'umor vitreo colma il volume tra il cristallino ed il fondo dell'occhio, una sostanza 
gelatinosa che ne determina la forma e la consistenza, permette il passaggio della luce verso la 
retina e lo protegge attutendo gli urti. 

L'occhio svolge quindi la funzione di lente capace di mettere a fuoco con straordinaria precisione, 
ma il trasduttore segnale-sistema nervoso si trova nella membrana retinica situata sulla superficie 
posteriore. 

La retina è suddivisa in molteplici strati (fig.1.2). Dopo una sottile membrana di separazione 
con l'umor vitreo la luce incontra la fibra nervosa, le sinapsi ed i gangli che raccolgono i segnali 
dalla membrana per trasportarli fino al nervo ottico (innestato nel punto cieco dell'occhio). Più 
in profonditàà vi sono neuroni adibiti alla integrazione e regolazione degli stimoli provenienti dal 
penultimo strato, che ospita i trasduttori veri e propri. L'mner nuclear layer (INL) contiene i 
nuclei di cellule orizzontali, amacrine e bipolari, che elaborano e trasportano il segnale dai fotore- 
cettori verso i gangli. I plexiform layers interni ed esterni sono formati dai dendriti appartenenti 
ai nuclei interni, e dalle loro intricate connessioni con i livelli adiacenti. 
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Figura 1.2: A sinistra: Stratificazione della retina (luce dal basso). A destra: Disposizione controintuitiva degli strati 
retinici, con le cellule sensibili a contatto con l'epitelio retinico pigmentato. 



Il penultimo strato della retina è costituito dalle cellule fotosensibili vere e proprie, le quali 
traggono nutrimento dal sottostante Epitelio Retinico Pigmentato cui sono ancorate. Oltre a 
nutrire ed allontanare gli scarti attraverso la barriera ematoretinica, l'epitelio impedisce distorsioni 
diffusive grazie al suo colore scuro; rigenera i recettori fagocitando le cellule inattive e rendendo 
nuovamente disponibili i pigmenti. La coroide avvolge infine tutto l'occhio nutrendolo attraverso 
fitti vasi sanguigni. 

La disposizione degli strati retinici è controintuitiva: la luce deve attraversarli tutti prima di 
arrivare alla retina. Per usare un'analogia ingegneristica, l'evoluzione ha preferito disporre il 
sensore in fondo alla telecamera^ dopo tutti i circuiti ed i dispositivi di controllo. 

Questo si spiega [16] pensando al danneggiamento continuo cui sono sottoposti i recettori, a 
causa appunto della loro spiccata reattività alla luce. La cornea, il cristallino ed i primi strati 
della retina svolgono efficientemente il ruolo di filtri, rimuovendo le lunghezze d'onda pericolose 
contenute nella luce solare. Ma i fotorecettori hanno una vita breve e sviluppano in continuazione 
sottoprodotti tossici a causa dei cammini di decadimento imprevisiti che una reazione fotochimica 
in un sistema così complesso può percorrere. 

Grazie alla disposizione invertita le cellule sensibili si trovano a contatto con l'epitelio retinico 
pigmentato capace di neutralizzare rapidamente le tossine, metabolizzare efficacemente le cellule 
distrutte dalle radiazioni assorbendo quelle in eccesso, e rigenerare continuamente il labile tessuto 
fotorecettivo. 

Le cellule fotorecettrici si distinguono in due tipi: bastoncelli e coni. I primi (fig.1.3) rispondono 
anche a luci fievolissime, tanto che basta un solo fotone per l'attivazione e sono disposti prevalen- 
temente nelle regioni periferiche della retina in numero di circa 120 milioni. Offrono una visione 
monocromatica ed indistinta, poiché un singolo neurone è incaricato di sommare ed amplificare il 
segnale prodotto da molteplici bastoncelli. 
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I coni sono invece assai meno numerosi, intorno ai 6 milio- 
ni, e concentrati principalmente nelldi fovea centralis^ il punto 
della retina adibito alla visione cromatica distinta. La loro 
sensibilità è massima in tre porzioni dello spettro visibile, 
centrate su 565nm (rosso), 535nm (verde), 440nm (blu), ma 
decade rapidamente all'indebolirsi della luminosità: un cono 
può necessitare centinaia o migliaia di fotoni per attivarsi. 
Il segnale di ogni cono viene elaborato da un neurone dedi- 
cato, trasmettendo immagini più chiare e distinte al cervello 
rispetto a quelle ottenibili dai bastoncelli. 

I segnali generati dalle due cellule si escludono tra loro, 
poiché la normale luce diurna inattiva i delicati bastoncelli 
mentre quella crepuscolare e notturna non è sufficiente ad at- 
tivare i coni. A questo corrispondono due qualità di visione: 
quella notturna monocromatica e sfocata, ma sensibilisima 
a luci deboli; e quella diurna ricca di colori e dettagli minuti. 

Bastoncelli e coni sono essenzialmente neuroni, con un as- 
sone proiettato verso il nuclear layer esterno, ma dotati di 
un cilio cellulare di foggia caratteristica, tubolare o conica, 

da cui i rispettivi nomi. Questa propaggine contiene gli organelli che ospitano il fotorecettore 
invaginato in dischi di membrana cellulare, impilati uno sull'altro. 

I bastoncelli sono stati studiati e compresi per primi, e la nostra ricerca riguarda il pigmento 
attivo in queste cellule: la Rodopsina. Fu il premio Nobel George Wald, a partire dagli anni 
30' [7,53-57], a gettare le prime luci sul loro funzionamento. 

Lo stimolo luminoso chiude i canali Na'^ e Ca^^ e causa l'iperpolarizzazione della membrana 
cellulare, condizione in cui il bastoncello non è più in grado di rilasciare il neurotrasmettitore 
glutammato. Il cambiamento di potenziale viene tramesso attraverso la sinapsi alle cellule (soli- 
tamente bipolari) del nuclear layer interno, quindi percorre i gangli, le fibre nervose retiniche, il 
nervo ottico e raggiunge il cervello. 




Figura 1.3: Schema di una cellula baston- 
cello [5]. La luce proviene 
dall'alto. 



1.2 II ruolo della Rodopsina nella visione 



Le prime intuizioni sul pigmento contenuto negli occhi dei vertebrati risalgono al 1877, in base 
alle osservazioni di Boll [4] sul colore degli occhi delle rane. Egli fu il primo a suggerire che il 
colore rosso del fondo del bulbo non derivava dall'emoglobina, ma da un pigmento nuovo proprio 
della retina. Kùhne [26] ne raccolse immediatamente le osservazioni: coniò il nome Rodopsina 
dal greco rhodon (porpora) e opsis (vista) per indicare la nuova specie colorante; isolò la retina 
identificandone gli strati; riuscì ad estrarre il pigmento dissolvendo i bastoncelli con bile, e molto 
altro. 

Fu Kùhne a traghettare la ricerca fino al punto in cui la riprese Wald negli anni 30', che 



1.2 II ruolo della Rodopsina nella visione 

scoprì [53] il ruolo nel ciclo visivo del carotenoide retinale, ovvero la vitamina A, contribuendo 
poi a chiarire la biochimica della proteina Rodopsina (Rh) [58]. 

La Rodoposina (fig.1.4) è il recettore transmembrana a sette eliche accoppiato alla G-proteina 
eterotrimetica Trasducina. È formata dalla proteina Opsina legata al cromoforo retinale. 
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Figura 1.4: Rappresentazione della Rodopsina transmembrana, in rosso il cromoforo. A destra: legame imminico (base di 
Schiff protonata) con la lisina. [8] 



Ogni bastoncello contiene circa 25.000 molecole di Rh per /im^, e l'estrema sensibilità è merito 
della reazione a catena [28] che si innesca dall'interazione con anche un singolo fotone. 

Il gruppo prostetico 11-cis-retinale legato alla Rh ha lunghezza d'onda d'assorbimento massimo 
a 500nm; dopo l'eccitazione elettronica alimentata dal fotone, esso rilassa allo stato fondamentale 
per via non radiativa con una torsione verso la conformazione all-trans-retinale. 

L'isomerizzazione del retinale da cis a trans, o E/Z in alcune notazioni^, induce un cambia- 
mento conformazionale nell'intera proteina Opsina, non in grado di legare la forma trans (fig.L5). 
Il complesso transitorio in cui l'Opsina è ancora immobile, detto Batorodospina (Ba), si defor- 
ma rapidamente con l'espulsione del retinale. Nel processo l'Opsina assume tre conformazioni 
intermedie: Lumirodopsina, Metarodopsina I, Metarodopsina II per finire con l'Opsina libera. 
L'ultimo intermedio è il fattore di scambio del nucleotide guanina (GEF) per l'enzima accoppiato, 
Trasducina, che viene quindi attivato (fig.LG). 

La trasducina catalizza la trasformazione di centinaia di molecole di guanosina trifosfato (GTP) 
in messaggeri difosfato (GDP), determinando una prima amplificazione del segnale. Il GDP a 
sua volta rompe il legame tra l'enzima Fosfodiesterasi (PDE) con il suo inibitore: la potente 
attività enzimatica che ne consegue catalizza l'idrolisi di più di quattromila molecole di guanosina 
monofosfato ciclica (cGMP). 

Il cGMP è responsabile dell'apertura dei canali sodio e calcio che consentono la diffusione degli 
ioni verso l'interno della membrana cellulare: la fosfodiesterasi ne interrompe quindi la funzione 
idrolizzandolo. I canali ionici si chiudono e la membrana si polarizza immantinente: mentre le 
pompe sodio-potassio continuano ad espellere ioni alla stessa velocità, l'ingresso ne è ostacolato. 
Si accumula così la differenza di potenziale che Passone segnala ai neuroni dell'INL, suscitando lo 
stimolo nervoso. 



""^E, entgegen: opposto, geometria trans; Z, zusammen: insieme, geometria cis 
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Figura 1.5: isomerizzazione del retinale nella cavità opsinica, incapace di legare la forma ali-trans [40]. Conformazioni 
intermedie della Rodopsina prima di espellere completamente il retinale [21]. La Fotorodopsina è una forma 
calda di Batorodopsina, in cui il retinale è altamente distorto. La sua esistenza è tuttora discussa 




Figura 1.6: Cascata biochimica che sfocia nell'iperpolarizzazione della membrana cellulare del bastoncello. Step): Fotorea- 
zione da Rodopsina a Metarodopsina. 2) Attivazione della Proteina G. 3) Produzione del complesso inibitorio 
tra una subunutà a di G ed il GTP. 4) Inibizione della fosfodiesterasi PDF. 5) Il cGMP viene esaurito dalla 
PDF ed i canali sodio non possono più rimanere aperti. [28] 
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L'evoluzione ha selezionato per la funzione visiva la più rapida trasduzione proteica nota, iniziata 
da una delle cinetiche fotochimiche più veloci ed efficienti mai scoperte in natura. L'intera cascata 
biochimica richiede circa 500ms [27]: il mutamento conformazionale da Rh a Ba ha luogo nei 
primi sei picosecondi e, di questi, solo 200 femtosecondi sono impiegati nella fotoisomerizzazione 
E/Z [47]. 

Per comprendere l'istantaneità dell'evento fotochimico se ne rapporti la durata a circa tre minuti 
(lfs=ls): sulla stessa scala dei tempi un battito del nostro cuore dovrebbe durare circa 27 milioni 
di anni! Indagini sperimentali con tale risoluzione temporale divennero praticabili solo dalla metà 
degli anni 80' del secolo scorso grazie alla elevatissima precisione raggiunta dalla tecnologia laser. 
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Figura 1.7: Spettroscopia laser al femtosecondo con metodo pump-probe (pompa-sonda) su un getto molecolare. [64] 

Fu il premio Nobel (1999) Ahmed H. Zewail ad applicare intensi impulsi di fotoni coerenti emessi 
da un mode-locked dye laser allo studio dello stato di transizione della dissociazione ICN->I+CN 
. Nel 1986 Zewail utilizzò un laser tradizionale {pump^) per innescare la reazione fotochimica in 
un getto molecolare, ed un laser sonda {probe) costituito da un treno di impulsi a distanza di 
pochi femtosecondi per interrogare il sistema (fig.1.7). 

Le conclusioni di questi esperimenti pump-probe portarono ad un ribaltamento di paradigma: 
non solo l'ineffabile stato di transizione diventava finalmente osservabile, ma nei femtosecondi 
successivi al cambiamento di stato le molecole si muovevano in modo coerente (fig.1.8). 

La nuova metodologia fu chiamata spettroscopia dello stato di transizione nel femtosecondo 
(FTS). 

Nel 1988 Zewail condusse l'esperimento destinato ad assumere lo status di "paradigma" della 
femtochimica: la FTS fu applicata alla dissociazione dello ioduro di sodio, un sistema con due 
potenziali (ionico e covalente), una coordinata nucleare e due stati elettronici. Le osservazioni 
furono accolte con grande eccitazione poiché riuscirono effettivamente a seguire la risonanza tra 
legame ionico e covalente, e l'avvicinarsi e allontanarsi periodico di sodio e iodio sullo stato 
eccitato. 

Il segnale prodotto dal laser sonda non poteva essere frutto del caos molecolare, giacché aveva 



^Laser "pompa": eleva le molecole dallo stato fondamentale a quello eccitato, come una pompa con l'acqua 
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Figura 1.8: A sinistra: stato eccitato e fondamentale per la molecola NaI. A destra: Oscillazioni coerenti della distanza 
internucleare Na-I sulla superfìcie di potenziale dello stato eccitato. [64] 



una precisa periodicità perfettamente distinguibile dal rumore di fondo. Questo portò ad una 
conclusione fondamentale: 

• L'eccitazione laser istantanea porta tutte le molecole coinvolte in conformazioni geometriche 
dalla distribuzione spaziale simile a quella dello stato fondamentale, in accordo con il prin- 
cipio di Franck Condon. Se nello stato fondamentale le molecole sono fortemente localizzate 
attorno alle posizioni di equilibrio, l'impulso istantaneo le pomperàà allo stato eccitato senza 
variarne le geometrie. Questo fenomeno è detto coerenza e permette di trattare la dina- 
mica di milioni di molecole mediante la traiettoria classica di una singola molecola, almeno 
finché perturbazioni Inter ed intra molecolari non inducono decoerenza. 



Nel 1990 Schoenlein, Shank et al. [46] costruirono un laser blu-verde capace di emettere impulsi 
ogni lOfs. Combinarono un laser controlled-pulse-mode-locked capace di emettere impulsi da 50fs 
con due amplificatori a colorante {dye laser) in serie, a loro volta alimentati da laser ad eccimeri di 
Xe-F. Dal raggio prodotto dal primo amplificatore veniva generata una luce con spettro continuo 
focalizzandolo su un getto di glicole etilenico; da questa sorgente si filtrava la regione d'interesse 
per deviarla nel secondo amplificatore. Da 50fs gli impulsi erano poi compressi a lOfs con tecniche 
ottiche. 

In questo modo Shoenlein poteva generare impulsi tra i 450nm ed i 580nm: l'anno successivo [47] 
li adoperò per dipanare l'evoluzione conformazionale della Rodopsina, con Xmax appunto nella 
regione blu-verde dello spettro. 

Gli spettri di assorbimento (fig.1.9) risolti nel tempo a 500nm mostrano la formazione dello 
stato eccitato al tempo (picco negativo), e la sua scomparsa in 125fs con recupero parziale in 
250fs ed esaurimento completo in 8ps. Gli spettri a 550nm e 570nm riportano l'assorbimento 



1.3 Femtochiinica della Rodopsina 
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Figura 1.9: Sinistra: cambiamento di assorbimento (ST/T) a 4 lunghezze d'onda. Destra: spettri differenziali in 4 tempi. 
[47] 

dei fotoprodotti, con picco a 200fs e poi stabilizzazione. Dagli spettri differenziali si identificano 
invece lo stato eccitato (490nm-540nm) dopo 33fs, ed i fotoprodotti dopo 200fs (570nm). 

Una velocità così elevata suggerisce l'esistenza di uno stato vibrazionalmente coerente e di 
una transizione senza barriera prodotto-reagente, con stati fondamentale ed eccitato fortemente 
accoppiati. 

Recentemente [25] la spettroscopia Raman risolta nel femtosecondo (FSRS) ha permesso di 
approfondire la dinamica. Il campione viene illuminato da un impulso Raman ogni picosecondo, 
mentre il continuum con periodo di 25fs sonda il sistema tra un impulso e l'altro. 

Nel corso di Ips la Rh si trasforma in Ba (fig.1.10), con red-shifting della banda 154:8cm~^— > 
1535c772~-'^ (stretching etilenico C=C); modificazione della zona 1100— ISOOctti"-'^ (stretching C-C); 
ed infine una inaspettata evoluzione della regione 800 — 950c7tz~^, prodotta dal wagging fuori dal 
piano dell'idrogeno (HOOP). 

Quest'ultima osservazione ha spinto alla formulazione di un modello teorico in cui la Rh si 
muove su una superficie di energia potenziale lungo gli HOOP per i primi 50fs, quindi l'energia si 
distribuisce anche lungo la torsione ed in Ips giunge alla Ba. 



1.3.1 Previsione della dinamica molecolare con metodo misto QM/MM 

Una previsione accurata dei primi istanti della fotoisomerizzazione della Rodopsina può essere 
ottenuta per via teorica con calcoli di dinamica molecolare. 

L'evoluzione cinetica di grandi molecole, come le proteine, è solitamente simulata risolvendo 
le equazioni di Newton per il moto in un campo di forza parametrizzato secondo i metodi della 
meccanica molecolare (MM). Per le reazioni fotochimiche la teoria MM non è applicabile: la 
parametrizzazione del campo di forza è infatti valida per un solo stato elettronico, solitamente 
quello fondamentale, mentre una fotoreazione implica almeno una transizione tra due stati. 
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Figura 1.10: A sinistra: Spettro Raman a tempi diversi. A destra: Modello dinamico di Kukura [25] 

I metodi di dinamica quantistica, basati sui principi della meccanica quantistica (QM), trattano 
in modo completo l'evoluzione della reazione sugli stati elettronici, ma il costo computazionale 
diventa proibitivo per molecole con più di qualche decina di atomi. 

Finora l'approccio a questo tipo di problema è stato studiare separatamente con metodo QM 
il cromoforo Retinale, o specie biomimetiche ancora più semplici. 

Martinez et al. [3] hanno spiegato la selettività con cui il retinale isomerizza esclusivamente lungo 
uno specifico doppio legame C=C tra i 5 possibili analizzando le superfici di energia potenziale 
dello stato eccitato e fondamentale della molecola ottenute dalla QM. 

Olivucci et al [13] hanno ricostruito il cammino di energia minima per la fotoisomerizzazione 
del Retinale con calcoli QM ab initio^. Hanno così potuto dimostrare che la reazione procede per 
un certo tempo di induzione lungo una coordinata collettiva di stretching, e solo successivamente 
isomerizza lungo quella torsionale (fig.1.11) - contrariamente a quanto ipotizzato in base agli 
esperimenti di spettroscopia laser esposti nel paragrafo precedente. Le conclusioni computazionali 
sono ribadite anche in [36,38,43,44] e confermate nella presente Tesi. 

Le pubblicazioni sullo studio del retinale e affini con metodi QM sono molteplici, ma chia- 
ramente possono dare una indicazione solo qualitativa sulla dinamica della proteina intera. In 
collaborazione con Manuel Domer e Daniele Passerone abbiamo scelto di sviluppare un meto- 
do misto QM/MM per simulare la dinamica della Rodopsina. Il campo di forze generalizzato 



"^A partire dai principi primi della fìsica e della meccanica quantistica, senza introdurre parametrizzazioni ottenute 
sperimentalmente a posteriori 
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1.4 Motori molecolari e switch fotochimici 
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Figura 1.11: Cammini di energia minima sullo stato eccitato per il retinale [13] 

GAFF [59] (MM) è stato parametrizzato in modo da includere il potenziale su stato fondamen- 
tale e primo eccitato del CFP (QM). Sono così stati condotti i primi calcoli di dinamica mista 
QM/MM e la consistenza dei risultati è stata valutata nel lavoro di tesi di Domer [6]. 

1.4 Motori molecolari e switch fotochimici 

L'elevata velocità della reazione fin qui descritta (<200fs) permette un'ottima resa quantica di 
0.67 [35], ovvero più di metà dei fotoni incidenti risulta reattivo con il 60% dell'energia trasferita 
al primo intermedio termodinamicamente stabile (Ba). 

Sistemi biomimetici del retinale sono assai appetibili per la costruzione di interruttori molecolari 
(switch), ovvero specie bistabili intercambiabili reversibilmente e non spontaneamente tra due 
o più stati [60]. L'attuazione risponde a stimoli, talvolta combinati, come la variazione di pH, 
concentrazione, temperatura; l'applicazione di un potenziale elettrochimico; l'esposizione alla luce. 

I primi impieghi tecnologici degli interruttori molecolari riguardano gli indicatori di pH e di 
concentrazione (sensori) a fini analitici; a partire dagli anni 70' lo sviluppo delle nanotecnologie 
ha impresso un forte impulso nella ricerca di nuovi switch, con l'obiettivo di realizzare motori 
molecolari e dispositivi supramolecolari capaci di svolgere funzioni complesse a livello nanometrico. 

II cambiamento di stato dello switch può alimentare un movimento articolato in un intorno mo- 
lecolare opportunamente ingegnerizzato (fig.1.12): in una affascinante rassegna introduttiva [1] 
sono illustrati molteplici esempi di "semplici" motori supramolecolari. Tra gli obiettivi a lungo ter- 
mine, la totalizzazione di lavoro macroscopico utile organizzando miliardi di movimenti collettivi 
nanometrici (ad esempio su superfici) e l'elaborazione dati parallela in computer molecolari [17,34]. 

Tra i possibili vettori utilizzabili per impartire energia e far scattare questi interruttori, la luce 
è il più promettente: 

• Il trasferimento energetico fotone -^ molecola è istantaneo; è richiesto solo il tempo per la 
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r igura 1.12: A sinistra: strutture complesse in grado di compiere semplici movimenti quando viene attivato uno switch 
presente sulla catena o sugli anelli [1]. A destra: braccio molecolare con estensione variabile basato su un 
diarilidene ed implementazione in un controllo di giunzione a livello nanometrico [9]. 

riorganizzazione conformazionale. L'immediatezza è un notevole vantaggio rispetto a switch 
termici o chimici che necessitano deha molto più lenta diffusione di calore o specie chimiche. 

• Nessuno scarto viene prodotto e nessun reagente deve essere ripristinato ad ogni cambia- 
mento di stato. 

• La reazione avviene contemporaneamente in tutto il volume in cui si è disposto lo switch 
(quelli elettro/chimici diffondono dall'elettrodo o dal punto di immissione dei reagenti) 

Alla classe dei fotoswitch appartengono azobenzeni, diarileteni/diarilideni, di-tienil-eteni, ani- 
dridi di-alchilidene-succiniche, stilbeni, quinoni fenossiantr acenici, spiropirani. Gli switch foto- 
chimici sintetici basati sull'azobenzene sono tra i più studiati in letteratura [41], per controllare 
la formazione di complessi con ioni; ripiegamento, poprietà elettroniche e catalitiche di catene 
peptidiche. 

Lo studio approfondito delle superfici di energia potenziale di uno switch permette sia la pro- 
gettazione di specie migliori, ad esempio rimuovendo barriere o inserendo vincoli sterici e repulsivi 
in modo da ottenere il movimento desiderato. Una possibilità molto innovativa è quella di dirigere 
le reazioni fotochimiche modulando accuratamente gli impulsi laser in modo da accompagnare le 
molecole lungo certi percorsi di reazione in modo assai specifico [29,45]. 

1.5 Obiettivi della Tesi 

I traguardi posti innanzi a questo lavoro di Tesi possono essere riassunti in tre macrobiettivi 
corrispondenti a tre periodi di ricerca ed articolati nei passaggi descritti di seguito e schematizzati 
in figura 1.13. 

1.5.1 Campionamento delle superfici di energia potenziale 

II primissimo passaggio per studiare il nostro sistema col metodo ab initio CASSCF è stata la 
definizione delle coordinate di reazione, e la preparazione di una OJtattice Z opportuna per il 
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programma di chimica quantistica Qaussian [11]. È stata necessaria una riflessione profonda 
sulla natura della reazione che saremo andati a descrivere con quella matrice Z (sez. 2.1.4). Ciò 
ha condotto alla scrittura di un programma dedicato alla manipolazione della OJlairice Z {GeoM<n 
, app. B) capace di muovere la molecola in modo chimicamente sensato. 

È seguita la mappatura delle superfici di energia potenziale (PES) su cui si muove il CPP negli 
stati fondamentale ed eccitato: stabilita una griglia di punti nelle coordinate di interesse si è 
svolto un calcolo CASSCF su ciascuno di essi. 

Le particolari esigenze del calcolo hanno impedito di avvalerci degli algoritmi di scansione già 
disponibili in Qaussian : abbiamo allora ideato il programma ReQuieM (app. A), in grado 
di percorrere la griglia in modo appropriato e di suddividere il carico computazionale tra diversi 
calcolatori e cluster con notevole accelerazione della scansione. 

Le coordinate importanti, identificate in letteratura e da considerazioni chimiche preliminari, 
apparivano essere quattro: una torsione, due piramidalizzazioni ed una coordinata di stretching 
lungo i doppi legami. Prima di affrontare una scansione quadridimensionale con tutti i pro- 
blemi computazionali e interpretativi che ne sarebbero conseguiti, abbiamo esplorato lo spazio 
tridimensionale di torsione e due piramidalizzazioni. 

Questi primi risultati hanno indicato che una delle due piramidalizzazioni ha effetti trascurabili 
sulle PES, ed hanno fatto emergere una chiara correlazione tra la torsione e la rimanente piramida- 
lizzazione. Analizzando gli orbitali dello spazio attivo in alcuni punti dalle coordinate caratteristi- 
che è stato possibile formulare una ipotesi sull'origine dell'interazione torsione-piramidalizzazione. 
L'ipotesi è stata convalidata mediante un modello ridotto con orbitali di Hùckel. 

Scartata dunque una coordinata di piramidalizzazione, abbiamo definito quella di stretching 
unificato e realizzato una seconda griglia tridimensionale con-^LCAS [22] per ottenere, oltre alle 
energie, anche i momenti dipolari punto per punto e le cariche. 

1.5.2 Ricerca di una forma funzionale fisicamente significante per le PES 

Per introdurre l'informazione fisica proveniente da questi calcoli in un programma di dinamica o 
meccanica molecolare, evitandone però il costo computazionale, non basta una semplice interpo- 
lazione delle PES, poiché questi programmi debbono conoscere anche l'accoppiamento adiabatico 
tra gli stati elettronici e la forma delle superfici diabatiche. 

Gli stati elettronici adiabatici, dipendenti in modo discontinuo dalle coordinate nucleari intorno 
all'intersezione conica, sono stati proiettati su una base diabatica, ottenendo stati diabatici che 
variano in modo smooth in funzione della geometria. Una prima diabatizzazione con il metodo 
di Stock [18] ha prodotto risultati numericamente accettabili per le energie, ma non compatibili 
con l'andamento delle cariche e dei momenti dipolari lungo le tre coordinate. Abbiamo quindi 
scelto l'operatore monoelettronico dipolo elettrico quale grandezza adatta per costruire funzioni 
elettroniche diabatiche, calcolandone il valore in ogni punto della griglia cou-Wlcas 7. 

Abbiamo quindi diabatizzato le superfici adiabatiche partendo dai momenti dipolari, ottenendo 
due superifici diabatiche sostanzialmente differenti dalle precedenti. Queste sono state interpolate 
con una forma funzionale che ne assicura la corretta periodicità ed un comportamento accettabile 
entro la regione accessibile dal sistema. 
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Il fit delle superfici diabatiche è stato implementato in Fortran per l'integrazione con il campo 
di forze di meccanica molecolare sviluppato da Manuel Domer e Daniele Passerone per l'intera 
proteina, Rodopsina. 

1.5.3 Calcoli di dinamica semiclassica 

Una volta definite le PES adiabatiche e diabatiche diventa possibile simulare la dinamica classica, 
semiclassica e quantistica del CPP nelle prime centinaia di femtosecondi dopo lo stato eccitato. 

La dinamica semiclassica è stata calcolata implementando le forme funzionali per il potenziale 
sui due stati nel programma TSH , scritto da F. Santoro. Un pacchetto d'onda di 10^ traiettorie è 
stato fatto evolvere per 2 picosecondi sulle superifci adiabatiche integrando le equazioni di Newton; 
il surface hopping, cioè la transizione tra Si ed So è valutata dall'algoritmo fewest switches di 
Tully. 

Questi primi calcoli andranno in futuro integrati con la dinamica quantistica per definire un 
protocollo di validazione per metodi di surface hopping semplificati, applicabili alla dinamica di 
tutta la proteina. 
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Figura 1.13: Diagramma del lavoro di Tesi 
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2 Campionamento delle superfici di energia 
potenziale 



Nell'ambito dell'approssimazione di Born-Openheimer (BO) la costruzione della superficie di ener- 
gia potenziale procede fissando la geometria molecolare, espressa in coordinate^ a conformazioni 
diverse, punti, facenti parte dello spazio geometrico che si ritiene possa essere esplorato dalla mo- 
lecola durante il processo di fotoisomerizzazione. Si ottiene così una griglia di punti per ciascuno 
dei quali si calcola poi l'energia degli stati elettronici di interesse. 

La pur semplice molecola del 4-(cyclo-pent-2-enyldiene)-3,4-dihydro-2H-pyrrolinium+ (fig.2.1) 
possiede 3 * 22 — 6 = 60 gradi di libertà: una ipersuperficie di energia potenziale in un tal numero 
di dimensioni porrebbe difficoltà interpretative enormi, ma prima ancora, non potrebbe essere 
calcolata in tempi ragionevoli con le tecnologie attuali. Fortunatamente non è necessario studiare 
tutto lo spazio geometrico: la griglia includerà solo i gradi di libertà maggiormente coinvolti dalla 
fotoisomerizzazione, ovvero le nostre coordinate di reazione (di seguito chiamate Crea)- 
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Figura 2.1: La nostra numerazione del CPP (idrogeno e dummies esclusi). Conformazione trans. Gli atomi 3,2 sono detti 
interanello o giunzione. 



2.1 Le coordinate di reazione 

Alcune informazioni basilari sulle coordinate significative per la fotoisomerizzazione del CPP 
possiamo dedurle dal tipo di reazione in oggetto. Ad esempio, possiamo da subito osservare 
che l'angolo diedro tra i piani formati dai due anelli dovrà percorrere 180° per portare la molecola 
dalla forma Trans a quella Cis. 

Poiché l'asse torsionale è un doppio legame, questi dovrà rompersi e quindi allungarsi per 
consentire la torsione tra gli anelli. Il legame 3=2 fa parte di un sistema tt coniugato, quindi tutta 
la catena deve riorganzizzarsi per ridistribuire i due elettroni in eccesso, accorciando i legami 
singoli 7-3, 2-6 ed estendendo 11=7 e 6=10. 

A conferma di questa prima intuizione chimica, è stato proposto da diversi autori [3,18,25,36,42] 
su basi sperimentali e teoriche che la dinamica successiva alla fotoeccitazione del retinale e affini 
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2.1 Le coordinate di reazione 

(di cui CPP è modello) proceda prima lungo gli stretching coniugati, e solo in seguito lungo la 
torsione che conduce all'intersezione conica. 

Il sistema tt possiede altri due importanti gradi di libertà: le piramidalizzazioni degli orbitali 
sp2 degli atomi di giunzione 3 e 2. L'importanza della piramidalizzazione nella fotoisomerizzazio- 
ne di sistemi tt delocalizzati è stata evidenziata in [2,24,36]. Già nel 2000 Ben-Nun e Martinez [2] 
avevano calcolato le superfici di energia potenziale per l'isomerizzazione dell'etilene, attribuendo 
alla piramidalizzazione dei due atomi di carbonio un ruolo fondamentale per il raggiungimento 
dell'intersezione conica. Estendendo il modello nel 2002 Olivucci [36] ha svolto calcoli ab initio 
sul sistema tt isolato, il catione Z-penta-3,5-dieniminio, osservando che un segmento del cammino 
di fotoisomerizzazione presentava un accoppiamento tra torsione e piramidalizzazione dei due C 
sull'asse rotazionale. In [24] è stata poi studiata sperimentalmente e teoricamente la fotoisomeriz- 
zazione del 4-Me(C6H9)CHF, evidenziando una notevole correlazione tra la coordinata torsionale 
lungo la quale avviene la reazione, e quella di piramidalizzazione. 
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Figura 2.2: Sinistra: asse torsionale del Z-penta-3,5-dieniminio [36] / Destra: Coordinata di piramidalizzazione per il 
4-Me(C6H9)CHF [24] 

In base a queste informazioni definiamo un sistema di quattro coordinate di reazione. Crea- 

• Lo stretching (Z) esprime la variazione delle distanze di legame degli atomi del sistema 

TV 

• La torsione (X) è l'angolo di isomerizzazione 

• Le piramidalizzazioni di pentene (P2) e pirrolinio (P3), indicano la deformazione della 
geometria planare sp2 degli atomi 3 e 2. 

Per poter svolgere simulazioni di dinamica determiniamo la massa o il momento d'inerzia as- 
sociati a queste coordinate, assumendo in prima istanza che essa non dipenda dalle coordinate e 
che i termini di accoppiamento nell'Hamiltoniano cinetico siano trascurabili. 

Vale la pena di notare che queste approssimazioni si applicano solo alla dinamica esplorativa 
sul nostro modello a dimensionalità ridotta; infatti le traiettorie di dinamica molecolare QM/MM 
della proteina condotte dai nostri collaboratori [6] sono calcolate in coordinate cartesiane senza 
dunque introdurre alcuna approssimazione sull'Hamiltoniano cinetico. 
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2 Campionamento delle superfìci di energia potenziale 

2.1.1 Definizione dello Stretching, Z 

Da una precedente pubblicazione del nostro gruppo di ricerca [42] sono note alcune strutture 
rilevanti per la fotoisomerizzazione del CPP: la geometria di minimo energetico sullo stato fonda- 
mentale (trans), l'Intersezione Conica (CI), la geometria cis, e due stati di transizione (TS1,2) o 
minimi locali sullo stato fondamentale. 

Confrontando i due punti più significativi per la reazione in oggetto, cioè trans e CI, troviamo 
conferma delle ipotesi riguardo lo stretching: il legame 3=2 (B2) si allunga da 1.356 a 1.446; 
simile sorte spetta agli altri doppi legami, mentre i legami singoli si accorciano da 1.44 ad 1.39 
(dati in tab.2.1). 

11 ^15^7 5 ^ 



N ^ \B6 



2 



/ B2 \B5 
12 ^^ / \ ^^ 10 

Figura 2.3: Nomenclatura dei legami sul sistema n 

Chiamiamo le 5 distanze interatomiche Cz che legano il sistema tt come mostrato in figura 
2.1.1, e definiamo una coordinata di stretching alternato Z tale che: 

• Per Z=0, tutte le distanze assumono i valori della geometria (C) trans: 

• Per Z=l, assumono i valori della CI: 

Cd = {BWci. BQci. B2ci, Bbci, B9ci} 

• Per 0<Z<1 le distanze variano linearmente con Z tra il valore trans e quello CI; per Z>1 e 
Z<0 si estrapola lo stesso andamento. 

Introduciamo un vettore k dato dalla differenza tra Cd e Ctrans'- 

k = Cci - Ctrans ={BlQci " Bl{)trans, BQci " BQfrans, 

B2ci - B2transi B5ci - Bòfrans, B9ci - Bdfrans} (2.1) 

La coordinata di reazione Z definisce una geometria C{Z) secondo l'equazione della retta: 

C{Z) = Ctrans + Z^k (2.2) 

Il valore di Z può essere ricavato da una geometria appartenente alla retta C{Z) tramite: 

^ ^ jCjZ) - Ctrans) ' ^ /^ o\ 
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2.1 Le coordinate di reazione 

L'eq.2.3 perde di significato per qualsiasi geometria non appartenga a C{Z). 

In tabella 2.1 sono riportate le distanze di legame per le geometrie Trans e CI con le loro 
differenze, ovvero gli elementi k^ deireq.2.1. Il segno positivo di k^ per i legami doppi (B10,B2,B5) 
e negativo per quelli singoli (B6,B5) indica lo scambio di ordine di legame; la variazione di distanza 
più consistente (più del doppio della media degli altri) è sul legame interanello 3^2 attorno al 
quale il CFP isomerizza. 



Var 


Trans 


CI 


CI-Trans 


BIO 


1.28756 


1.32287 


0.03531 


B6 


1.43415 


1.39227 


-0.04188 


B2 


1.35570 


1.44601 


0.09031 


B5 


1.45774 


1.39907 


-0.05866 


B9 


1.34658 


1.37626 


0.02969 



Tabella 2.1: Distanze di legame sul sistema n e coefficient 



1 cz 



Massa associata allo stretching Z 

Lo stretching Z è una coordinata interna: per ricavare la massa associata al moto lungo Z abbiamo 
applicato il metodo esposto in "Molecular Vibrations" ( [61], cap. 4), dove sono trattate vibrazioni 
molecolari in coordinate interne. 

Consideriamo una distanza di legame s^^k G (B10,B6,B2,B5,B9) ed esprimiamola in funzione 
delle coordinate cartesiane di spostamento Xi degli atomi della molecola: 



s^ = Y^ B^Xi + 9 X] ^iJ^i^J + termini superiori 



(2.4) 



*,j 



I coefficienti B^^ B^ ■ sono determinati dalla geometria molecolare e stabiliscono come la lunghezza 



*,j 



di legame varia in funzione delle coordinate cartesiane. Se assumiamo variazioni di posizione molto 
piccole possiamo trascurare i termini di ordine superiore e scrivere req.2.4 in forma matriciale: 



s = 



(2.5) 



Includendo in s le sei coordinate mancanti per descrivere traslazione e rotazione, req.2.5 diventa 
invertibile: 



(2.6) 



Ora possiamo scrivere l'energia cinetica in coordinate interne: 



2r = x^Mx = 5^B-^tMB-^5 = 5^G"^5 



'^siJq 



k^Vk 



= mk 



(2.7) 
(2.8) 
(2.9) 



dove la matrice diagonale M contiene le masse nik dei 6 atomi e G è la matrice energia cinetica 
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2 Campionamento delle superfìci di energia potenziale 

in coordinate interne. Dalla matrice G si ricavano le masse associate agli stretching, ma prima 
occorre definire la matrice dei coefficienti B e quella delle masse M. 

Assumiamo che lo stretching avvenga su un piano, e impieghiamo così solo due dimensioni (x,y) 
per definire Xi- Trascuriamo inoltre gli accoppiamenti con tutte le altre coordinate interne. La 
matrice B ha sei righe: una per ogni lunghezza di legame r; e dodici colonne: due (x e y) per ogni 
atomo. I coefficienti Bij per ciascuna distanza di legame sono calcolati dalla posizione dei due 
atomi a,b tra cui è definita: 



B 



^k^Jxa 



rab 



B, 



B 



^k,3x 



rab 



^ab — Xa — X]j 
^ab — ^a — Xf) 



B, 



^kJya 



kJy^ 



rab 



Kb 

rab 



(2.10) 



dove ik è la riga corrispondente ad una lunghezza di legame 5^, jx^ è l'indice di colonna associato 
alla coordinata x dell'atomo a, jx^ per l'atomo b. 
La matrice B nell'insieme risulta: 



/ 



^1,2 
^1,2 








-^1,2 
^1,2 



^1,2 
■^1,2 


^1,2 
■^1,2 



























^2,3 


^2,3 


Ai 


3 


^'2 


3 

















^2,3 


^2, 


3 


^2, 


3 








Af 


4 


^"s 


4 


Af,4 
^^3,4 


A|,4 
^3,4 











^3, 


4 


rs, 


4 


















A|,5 
^^4,5 


^4,5 


A|,5 
^^4,5 


■^4,5 



























Af,6 
^5,6 


^5,6 


Af,6 
^5,6 



\ 



(2.11) 



La matrice M si costruisce secondo lo stesso schema. 

Calcoliamo G da BM~"'^B^. Considerando alcuni vincoli geometrici: 



(A-2)' + (AI2)' = rl2 

(Af,^)' + (AI,6)^ = ri, 

AÌ,3Af,4 + A|,3Al,4 27r 

— ■ ■ ■ = cos 



^3,4^2,3 

A a: A X I /\y A 2/ 
^5,6^4,5 ^ ^5fi^4 



4,5 



3 

27r 



cos 



(Ai,3)' + (Ala)' = rh 

i^h? + (Als)' = <5 

A-2Ai,3 + A|2A|3 2^ 
= cos — 

n,2^2,3 3 

^4,5^3,4 + ^4,5^3,4 ^TT 

— ■ ■ ■ = cos — 

^4,5^3,4 3 



^5,6^4,5 



si riduce a: 



G = 



/ i^c + m 


Mc/2 











i"c/2 


2|Uc 


-Aic/2 











-/^c/2 


2mc 


l^c/2 











/«c/2 


2i^c 


/^c/2 


V 








t^c/2 


2aìc 



(2.12) 
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2.1 Le coordinate di reazione 

dove /le = '^c^' I^N ~ ^iv^' ^^^ 1^ massa di carbonio e azoto in unità di massa atomiche a.m.u. 
Definiamo ora una base completa data da Z e 4 coordinate ortogonali Zg^s = 2,3,4,5. In 
tal modo abbiamo definito uno spazio in cinque dimensioni di cui gli Z^, dopo normalizzazione, 
costituiscono gli assi. Zi è il vettore delle differenze definito in tab.2.1, mentre gli Zg ortogo- 
nali sono stati calcolati da M. Doemer in [6] tramite una procedura Gram-Schmidt e risultano, 
normalizzati: 



Z = (0.281399, -0.337679, 0.723598, -0.474359, 0.241199) 
Z2 = (0.718875,0.694218,0.028819,-0.0188925,0.00960634) 
Z3 = (-0.0555096,0.0555096,0.588606,0.803436,-0.0432542) 
Z4 = (0.297245,-0.297245,-0.313073,0.313073,0.791999) 
Z5 = (-0.559104,0.559104,0.176362,-0.176362,0.559104) 

z = (Zi, Z2, Z3, Z4, z^y 



(2.13) 



Trasformiamo ora G nelle nostre coordinate Z, e sostituiamo le masse numeriche: 



& = Z^GZ = 

/ 0.140029 
-0.0264271 
0.029745 
-0.000986697 
V -0.0197066 



-0.0264271 

0.200375 

-0.0157034 

0.00668786 

-0.00127491 



0.029745 
-0.0157034 

0.200163 
0.0327323 
0.00353748 



-0.000986697 

0.00668786 

0.0327323 

0.162992 

0.0313799 



-0.0197066 \ 
-0.00127491 
0.00353748 

0.0313799 

0.117869 / 



Come massa associata alla coordinata Z scegliamo l'inverso del primo elemento diangonale di G', 
trascurando così le altre coordinate ortogonali, che rimangono fisse nel nostro modello, nonché 
qualsiasi accoppiamento: 



1836.1515aM 



mz = l/G'i 1 = 7.14138amu, 7.14138amu ■ 

' ^'^ ' 1.007825amu 

2.1.2 Definizione delle Piramidalizzazioni, P2, P3 



= 13010.8aM 



(2.14) 



Le analisi dei già citati [2,24,36] suggeriscono di costruire le 
coordinate di reazione intorno agli atomi 2 e 3 per descrivere 
la deviazione dalla planarità dei legami formati dagli atomi in- 
teranello con i tre atomi vicinali, indice del progressivo cambio 
di ibridazione sp2 -^ sp3. 

Inoltre in [36] è stata verificata per il CPP l'esistenza di una 
linea di intersezione conica, intersection seam, lungo la quale 
la differenza tra l'energia dei due stati è nulla e che porta dal 
punto Frank-Condon alla CI ed al Cis distorcendo la planarità 
degli anelli attorno agli atomi 2,3. 

Per ottenere una coordinata perfettamente indipendente da 



^t-z>P 




Figura 2.4: Definizione dell'atomo fitti- 
zio e della coordinata P 
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2 Campionamento delle superfìci di energia potenziale 

quella di stretching e di torsione è stato necessario introdurre 

due atomi fittizi 1 e 4 (dummy atoms) associati a quelli giunzione (2 e 3) definiti in modo che i 

tre angoli (P^) dummy-giunzione-vicinale siano identici. 

Dato un atomo giunzione G= {xgjVGj ^g) e tre atomi vicinali A,B,C, si definiscono: 

= (0,0,0) G = OG À=GA B = GB C = GC F = GF 
G-A 



Ra = 
Rb = 
Rc = 



\G-A\ 
G-B 
\G-B\ 
G-C 



\G-C\ 
d = \G — F\ = 1 distanza arbitraria 

L'incognito atomo fittizio F associato a G dovrà soddisfare il seguente sistema di equazioni: 

Ra ■ Rp = cos Pg 
R-B • Rp = cos Pg 
Re ■Rf = cos Pg (2.15) 

Il sistema può essere risolto in forma matriciale introducendo la matrice S degli atomi vicinali 
e K dei termini noti, quindi moltiplicando l'inverso di S per K, come mostrato in eq.2.16. 



K= cosPn ¥ 






(2.16) 

Il sistema diventa irrisolvibile qualora A,B,C,G siano complanari; si applica allora banalmente: 

Ra X Rb 



F 



\Ra X Rb\ 



Le coordinate P hanno valore 90° per piramidalizzazioni nulle; <90° per distorsioni verso il 
dummy; >90° per distorsioni opposte (in fig.2.4 è mostrato un angolo P>90°). 

L'intersection seam presente nella molecola a 3N-6 gradi di libertà lungo le piramidalizzazio- 
ni, viene ben approssimato dal nostro modello con linee di quasi-degenerazione in prossimità 
dell'intersezione conica, vedi sez.2.4.3. 
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2.1 Le coordinate di reazione 

Momento d'inderzia associato alla piramidalizzazione P3 

Il momento d'inerzia associato alla coordinata di piramidalizzazione è il momento d'inerzia ridotto 
calcolato per una rotazione rigida attorno ad un asse perpendicolare a 3,2. L'asse è complana- 
re all'anello del pentene e passante per 2 per la coordinata P2; giace invece sullo stesso piano 
dell'anello pirrolinico e passa per 3 per la coordinata P3. 

I due momenti d'inerzia Z^^"", /^ ^"^ sono calcolati l'uno per l'anello cui appartiene n per cui passa 
l'asse di rotazione Axn] l'altro per l'anello per cui non passa l'asse: a, b G (pi=pirrolinio,pe=pentene). 
I momenti d'inerzia ridotti sono allora: 



pi 




pi 




tAx2 
^pi 




pe 




pe 


Ipe ^ — Z_^ ^i ^i,Ax2 ' 


tAx2 

Jpe 


= ^mi tIa^^ 


7-^X2 T^X2 
jP2 _ pe ^Pi 






tAx2 I tAx2 

J-pe -h Ipi 






tAx3 7^X3 
jP3 _ P^ pi 

tAx3 i tAx3 

J-pe -h Jpi 







(2.17) 

Solo il momento d'inerzia per l'asse P3 è stata calcolata per motivi che verranno spiegati in 
sez.2.4.5. 

7^3 = QQ^nau (2.18) 

2.1.3 Definizione della Torsione, X 

È la coordinata di isomerizzazione: all'intersezione conica (CI) i due anelli giacciono su piani 
perpendicolari, e la loro planarità è pressoché conservata - come in trans ed in cis dove sono 
inoltre complanari. Questo suggerisce che la coordinata di torsione dovrebbe essere in grado di 
muovere gli anelli intatti^ senza possibilità di introdurre o dipendere da eventuali distorsioni. 

Come mostrato in fig.2.5 X è definita come l'angolo diedro 1-2-3-4, ovvero l'angolo tra il piano 
in cui giacciono 1,2,3 e quello passante per 2,3,4. 

La coordinata di torsione X sarà dunque: 

• X=180° per la geometria trans, il minimo energetico allo stato fondamentale 

• X=0° in geometria cis 

• X^90° alla intersezione conica 

• Nel corso della isomerizzazione andrà dunque da 180° a 0° 

Momento d'inerzia associato alla torsione X 

Il momento d'inerzia ridotto associato alla torsione è calcolato per la rotazione dei due anelli 
attorno all'asse 2,3. Come nel caso della piramidalizzazione, definiamo due momenti d'inerzia 
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2 Campionamento delle superfìci di energia potenziale 




Figura 2.5: Proiezioni di Newman raffiguranti l'angolo di torsione. Da sinistra: geometria trans, CI, angolo X 

relativi ai due anelli separati, ma in questo caso rispetto ad un asse soltanto (2,3), quindi ricaviamo 
il momento d'inerzia ridotto: 



pi 



pe 



J-r)P. J-r) 



J-TÌP, \ J-TÌ 



pe ~r j-pi 



(2.19) 



Risulta per la torsione un momento d'inerzia ridotto di: 

I^ = 260.9ai^ 



(2.20) 



Anche per Ix, come già per Ipg, si è trascurata la dipendenza dalle altre coordinate di reazione 
per mantenere semplice il modello. 

2.1.4 Architettura della Orlatrice Z 

Il formato predefinito da Qaussian per esprimere la geometria molecolare con coordinate interne 
è la QJlatrice Z . Ciascuna riga descrive un atomo, al quale ci si riferisce con il numero di riga 
A^ (riferimento). Le posizioni relative sono individuate da variabili di tre tipi: distanze tra due 
atomi, angoli tra tre atomi e diedri tra quattro. Una riga è composta da una etichetta atomica 
nella prima colonna (C,N,H, etc) e fino a tre coppie (riferimento ad atomo precedente A^-^; nome 
della variabile) nelle colonne successive, come segue: 

1. Riga 1: contiene l'etichetta dell'origine Ai, e nessun'altra colonna. 

2. Riga 2: etichetta dell'atomo A2; riferimento all'origine Ai, variabile distanza tra Ai e A2. 

3. Riga 3: etichetta dell'atomo 3; riferimento a 1 o 2 A^, distanza 3-A5; riferimento all'atomo 
rimanente (Ac), variabile angolo tra gli atomi 3 - A^ - A^. 

4. Riga N: etichetta dell'atomo N; riferimento a qualsiasi a<N-l (A^), distanza N - A^; riferi- 
mento ad altro atomo (A5), angolo tra gli atomi N - A^ - A5; riferimento Ac, variabile angolo 
diedro tra N - A^ - A5 - Ac- 

Specificati tutti gli atomi con le loro variabili distanza, angolo, diedro, sono poi elencate tutte 
le variabili con i loro valori, che determinano di fatto la geometria della molecola. In tabella 
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A Et 


E 


istanza 


Ango 


lo 


Diedro 


1 X 












2 C 


1 


Bl 








3 C 


2 


B2 


1 


A2 




4X 


3 


B3 


2 


A3 


1 D3 


5 C 


2 


B4 


1 


A4 


3 D4 


eC 


2 


B5 


1 


A5 


3 D5 


7C 


3 


B6 


4 


A6 


2 D6 


sC 


3 


B7 


4 


A7 


2 D7 


9C 


5 


BB 


2 


A8 


6 D8 


10 e 


6 


B9 


2 


A9 


5 D9 


11 N 


7 


BIO 


3 


AIO 


8 DIO 


12 C 


8 


Bll 


3 


Ali 


7 Dll 


13 H 


5 


B12 


2 


A12 


6 D12 


14 H 


5 


B13 


2 


A13 


6 D13 


15 H 


9 


B14 


5 


A14 


2 D14 


16 H 


9 


B15 


5 


A15 


2 D15 


17 H 


10 


B16 


6 


A16 


2 D16 


18 H 


6 


B17 


2 


A17 


5 D17 


19 H 


7 


B18 


3 


A18 


8 D18 


20 H 


11 


B19 


7 


A19 


3 D19 


21 H 


12 


B20 


8 


A20 


3 D20 


22 H 


12 


B21 


8 


A21 


3 D21 


23 H 


8 


B22 


3 


A22 


7 D22 


24 H 


8 


B23 


3 


A23 


7 D23 




Figura 2.6: I primi 8 atomi definiti nella 
VJtaitict Z contengono la coordi- 
nata di torsione e le due di 
piramidalizzazione 



Tabella 2.2: Mainct Z per il CPP 

2.2 è riportata la OJlaitice Z che abbiamo scritto per il CPP (omettendo la l'assegnazione delle 
variabili). 

La matrice è stata progettata per separare al meglio le 4 coordinate di reazione. La torsione è 
indicata dal diedro D3; P3 è uno qualsiasi dei tre angoli equivalenti A3,A7,A6 e, similmente, P2 
è uno tra A2,A4,A5. Z è composta dalle distanze B10,B6,B2,B5,B9 come spiegato in 2. LI 

Ogni atomo nella matrice necessita di 3 punti di riferimento, che abbiamo scelto tra quelli 
utilizzati per definire le coordinate in modo da non introdurre distorsioni indesiderate negli anelli. 
Pertanto gli atomi degli anelli sono definiti rispetto a quelli di giunzione ed adiacenti: in tal modo 
torsione e piramidalizzazione muovono gli anelli intatti, assieme agli atomi di idrogeno (fig.2.L4). 

La 9Jtatrice Z così costruita è stata sottoposta ad alcune verifiche prima di utilizzarla nella 
costruzione di griglie di punti. 

2.1.5 Separazione torsione-piramidalizzazione 

L'introduzione degli atomi fittizi e le relative coordinate P3^2 e X sono state verificate su alcune 
strutture generate da Sinicropi et al. [42]. Al fine di descrivere la linea d'intersezione (intersection 
seam) erano stati determinati punti di degenerazione Si /So al variare dell'angolo diedro Ada=8- 
3-2-5, in modo da introdurre una torsione tra i due anelli. Le strutture erano poi state ottimizzate 
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2 Campionamento delle superfìci di energia potenziale 

per localizzare una struttura in cui Si e So fossero degeneri, mantenendo però congelato detto 
diedro Ada. 

Abbiamo tradotto quelle conformazioni nella nostra OJlaitice Z , calcolando le nuove coordinate 
Creai riportate in tabella 2.3. 



Ada° 


X° 


P3° 


P2° 


B2 A (Zb2) 


30° 


76.800976 


104.309053 


80.79383 


1.479194 (1.3674) 


50° 


83.072541 


100.282872 


83.37497 


1.466258 (1.2242) 


70° 


85.978797 


96.14375 


88.08703 


1.448354 (1.0259) 


90° 


89.300419 


90.371070 


90.736932 


1.446953 (1.0104) 


110° 


92.187290 


84.888504 


93.742263 


1.442436 (0.9604) 


130° 


95.286704 


80.077191 


97.473134 


1.465927 (1.2205) 


150° 


102.254183 


76.386980 


99.823974 


1.485567 (1.4380) 


CI: 87.431° 


89.536010 


89.871629 


90.989663 


1.44601 (1) 


IVans: 179.948° 


179.983774 


89.990184 


90.007949 


1.35570 (0) 



Tabella 2.3: Espressione nelle nuove coordinate Crea- Lo stretching non può essere definito, quindi è riportato il valore 
del legame interanello B2 e l'inverso di eq.2.2 

Si nota che le strutture lungo l'intersection seam sono caratterizzate da coordinate X, P2 e 
P3 con variazioni decisamente inferiori rispetto a quella dell'angolo Ada: nel calcolo con vincolo 
Ada=30°, cioè -60° rispetto alla CI, l'ottimizzazione porta la torsione a -14° e le piramidalizzazioni 
a P2+14° e P3-10°. Con il vincolo Ada=150° le differenze con le coordinate della CI di minima 
energia sono anche minori. 

L'ottimizzazione vincolata dell'intersezione conica, dunque, riarrangia la torsione e le pirami- 
dalizzazioni come da noi definite in modo da portarle verso valori prossimi alla CI di minima 
differenza energetica tra i due stati elettronici, ovvero X e P2,3=90°, pur conservando l'angolo 
Ada. 

Questa ridistribuzione della deformazione tra X e P2,3 testimonia l'efficace disaccoppiamento 
tra le tre coordinate da noi scelte, e conferma il senso fisico dei loro punti di equilibrio: possiamo 
supporre che quanto più la molecola si avvicina a X=90°,P2,3=90°, tanto più si approssima ad 
una intersezione conica. 

2.1.6 Correzione dei diedri esterni agli anelli 

L'esplorazione delle geometrie al variare di torsione e piramidalizzazione ha permesso di indivi- 
duare un difetto importante nella OJlairice Z . 

In figura 2.7 è riportata la struttura trans con una piramidalizzazione su 2, P2=110°. Si osservi 
la deformazione del legame 9-10, che dai 1.50548A a P2=90° passa a soli 1.14925A: ciò è dovuto 
alla definizione della OJlairice Z . 

Essa è concepita per esprimere le posizioni atomiche utilizzando solo le variabili indispensabili, 
e la distanza 9-10 è ridondante visto che sia l'atomo 9 che l'atomo 10 sono già definiti (9 rispetto 
a 5,2,6; 10 rispetto a 6,2,5). Tale distanza viene dunque calcolata conseguentemente alle posizioni 
di questi atomi, ed in particolare il problema riscontrato tra 9,10 può essere ricondotto agli atomi 
5 e 6 che compaiono nelle definizioni di tutti i C dell'anello. 
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Figura 2.7: Distorsione sul legame 9-10 provocata da P2 

L'atomo 5 è definito dalla distanza B4 rispetto al 2, dall'angolo A4(=P2) formato con 2 ed il 
dummy 1, ed infine dal diedro D4 tra il piano 5,2,1 ed il piano 1,2,3. Il 6 è definito con la distanza 
B5 sempre rispetto a 2, l'angolo A5=P2 ed il diedro D5 tra il piano 6,2,1 e ancora il piano 1,2,3. 

I diedri D4 e D5 sono detti angoli diedri esterni agli anelli^ o solo diedri esterni. Ciò che accade 
quando si stabilisce P2 7^ 90° è che gli atomi 5 e 6 si debbono trovare sempre su un piano avente 
angoli diedri D4 e D5 fissi rispetto al piano 1,2,3, come schematizzato in figura 2.8. 

S'immagini l'atomo 2 al centro della sfera, la quale rappresenta la condizione di distanza di 
5 e 6 dall'atomo 2 ponendo che sia uguale per i due atomi, ovvero B4=B5. Il piano verde è il 
riferimento 1,2,3 mentre in blu e rosso sono disegnati i piani su cui gli atomi 5 e 6 sono vincolati 
dai diedri D4 e D5. Gli atomi 5 e 6 possono allora trovarsi solo sulle curve di intersezione tra i 
piani rosso e blu e la sfera, a seconda dell'angolo di piramidalizzazione. Per P2=90° l'angolo 5-2-6 
è massimo, mentre è nullo per P2=0°. È come se 5,6 si trovassero all'estremità delle stecche di 
un ombrello: quando è aperto l'angolo tra le stecche è massimo per tendere il tessuto; man mano 
che si chiude le stecche si avvicinano fino a toccarsi, e l'ombrello ritrova posto nello zaino. 

La contrazione dell'anello con avvicinamento di 9-10 non solo introduce un accoppiamento tra 
piramidalizzazione e stretching, ma si ripercuote sull'energia provocando un aumento ripidissimo 
in funzione di P2. 

Tuttavia nel corso della nostra reazione è ragionevole pensare che la distanza 9-10 non venga 
affatto coinvolta poiché estranea al sistema tt. Il problema viene risolto imponendo D4=D5 e 
ricalcolandoli in modo da mantenere fisso l'angolo interno 5-2-6 indefinito in filatrice Z . 

Consideriamo un sistema cartesiano con origine O, un punto A sull'asse z e tre punti B90, C90 
e D90 posti sul piano xy, tutti a distanza costante 1 da O; in aggiunta B si trova sul piano yz. 
Definiamo l'angolo P come AOC = AOB = AOD = 90°, il diedro come A = COAD = BOAD, 
e l'angolo interno a = BOC . Quanto deve valere A affinché, posto un certo P, l'angolo a sia 
ancora uguale ad ap=9o° 



? 



Supponendo di mantenere fisso il punto A, il cambiamento di P porta B90, C90 e D90 fuori dal 
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2 Campionamento delle superfìci di energia potenziale 




Figura 2.8: A sinistra: rappresentazione geometrica delle tre variabili che determinano la posizione degli atomi 5 e 6. A 
destra: sistema cartesiano per il ricalcolo degli angoli diedri esterni 



piano xy, nei punti B,C e D in fig.2.8. Troviamo il punto B, cui imponiamo di rimanere sul piano 
yz: 

= (0,0,0) 

A =(0,0,1) 

xb = (sul piano yz) 

B ' A = cos P ^ zb = cos P (angolo P) 

x\ + y\ + Zq = 1 ^ Vb = sin P (norma 1) 



B = (0,sinP,cosP) 



La posizione del punto C: 



C • A = cos P ^ zc = cos P (angolo P) 
B ' C = cos a (angolo a) 



Xc 

(0,sinP,cosP) • I yc 
cosP 



COSQ^ 
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2.1 Le coordinate di reazione 



cos a = yc sin P + cos^ P 
cos a — cos^ P 

yc = ^-B 

\C\ = 1 (norma 1) 



I 9 ^ COS a — cos^ P 
^C = A/ sm^ P 



sinP 

Proiettiamo ora B, C e l'incognito D ortogonalmente sul piano xy. Dxy dovrà rispettare la con- 
dizione DxyOBxy = DxyOCxy^ polclié dcvc glaccrc sullo stesso piano AOD rispetto al quale sono 
definiti i diedri uguagliati. Ricaviamo l'angolo ai compreso tra BxyOCxy e posizioniamo T)xy sulla 
bisettrice dell'angolo ai: 



Bxy = {xb, Vb), Cxy = (xc, yc) 
ai = arccos 



^xy ' ^xy 

^D^y = sin(7r + ai/2) 
VD^y =cos(7r + Q^i/2) 

Per ottenere D applichiamo la medesima piramidalizzazione P di B e C anche a Dxy^ quindi 
aggiungiamo z = cosP e moltiplichiamo x ed y per sinP (è una rotazione del vettore T)xy nel 
piano T)xyO^ di P). 

I nuovi angoli A=D4=D5, adeguati per mantenere l'angolo interno a fisso pur applicando una 
piramidalizzazione P, si calcolano infine dal prodotto scalare delle normali ai piani DOA e BOA 
(o equivalentemente COA): 

xd = XD^^sinP 
Vd =XD,ySmP 
zd = cos P 

m = e X A 

n2 = D X A 

ni '712 
A = arccos - — - — r 



FilF2| 



Si è fatto riferimento finora all'anello pentenico, ma il ragionamento è applicato anche alle 
coordinate della OJtairice Z che descrivono il pirrolinio. 
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2 Campionamento delle superfìci di energia potenziale 




Figura 2.9: Geometria ottenuta ponendo P2 = 110° e X=209° e applicando la correzione degli angoli esterni. Gli atomi 
d'idrogeno 14-19 sono molto ravvicinati. 

2.1.7 Approssimazioni per piramidalizzazioni elevate 

Si consideri ora la geometria con correzione degli angoli esterni P2=110° e X=209° riportata in 
fig.2.9. Gli atomi di idrogeno 14-19 si trovano ad una distanza di soli 1.62688A. Per testare gli 
ingombri degli H legati a 5,6 abbiamo calcolato la distanza tra ciascuna coppia appartenente ad 
anelli diversi (13,14,18 per il pentene; 19,23,24 per il pirrolinio) in funzione della torsione e per 
due forti piramidalizzazioni, P2=70° e P2=110°; in fig.2.10 sono riportate in grafico. 



Distanze H-H perP2 = 110 





Distanze H-H per P2=70 





Torsione X° 



Figura 2.10: Distanze H-H in funzione della torsione, calcolate per due piramidalizzazioni estreme diverse (110° e 70°) 

A piramidalizzazione molto maggiore di 90° l'idrogeno 14 si appressa agli H dell'altro anello 
al di sotto di 1.7A; per P< <90° è invece l'atomo 13 ad avvicinarsi troppo (seppur leggermente 
meno). 

Fortunatamnete l'analisi delle PES mostra che l'energia lungo le piramidalizzazioni non per- 
mette alla molecola di acquistare accelerazioni verso strutture talmente distorte da rendere si- 
gnificativo l'ingombro sterico degli H. Pertanto le esplorazioni di dinamica condotte con il nostro 
modello non risentono di questo effetto. 
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2.2 Orbitali dello Spazio Attivo 

Con una OJlaitice Z corretta è possibile descrivere qualsiasi posizione relativa degli atomi che 
contiene, ma l'obiettivo di una matrice ben architettata è quello di esprimere le coordinate di 
reazione (da noi definite) utilizzando il minor numero di variabili possibile (definite nella matri- 
ce). Nel nostro caso, la matrice ideale richiederebbe di modificare solo 4 variabili. Tuttavia non è 
possibile: la coordinata di stretching è collettiva e quella di piramidalizzazione introduce deforma- 
zioni complesse, che è necessario correggere modificando altre variabili ancora. Grazie ai metodi 
implementati nel programma GeoMan (app.B) scritto appositamente in linguaggio Python e Sci- 
Py [20,39], è possibile generare una geometria a partire dai 4 valori di Crea- H programma si cura 
di muovere in modo appropriato tutti gli stretching coinvolti, la torsione e le piramidalizzazioni 
assieme agli angoli esterni, come spiegato nella presente sezione. 

2.2 Orbitali dello Spazio Attivo 

Il metodo CASSCF richiede la definizione degli orbitali dello spazio attivo attraverso il quale 
verranno esplorate le eccitazioni elettroniche. Il processo in oggetto è innescato da un'eccitazione 
So ^ Si, quindi è presumibile che siano coinvolti gli orbitali tt occupati a più alta energia e quelli 
virtuali ad energia più bassa. Il sistema tt ha sei elettroni localizzati nello stato fondamentale in 
tre orbitali tt. Per descrivere lo stato eccitato includiamo nello spazio attivo anche i tre orbitali tt 
virtuali a più bassa energia. 

Per stabilire quali tra gli orbitali occupati appartegono al sistema tt abbiamo disegnato l'isosu- 
perficie di densità elettronica degli ultimi 5 orbitali occupati e dei primi 5 liberi, considerando la 
geometria Trans come riferimento (presentati in figura 2.11 ). 

Tra gli orbitali analizzati, risultano avere simmetria tt i 35(HOMO-l), 36(HOMO), 37(LUMO), 
38(LUMO+l). Il 34, naturale scelta per ordine energetico, deve essere scartato poiché ha quasi 
esclusivamente carattere a. L'orbitale ad energia più bassa ancora presenta lo stesso problema: 
la scelta ricade dunque sul 32(HOMO-4) in sostituzione del 34. 

Un ragionamento analogo si applica all'orbitale 39, il quale concentra la maggior parte della 
densità elettronica sull'azoto ed è anch'esso, per il resto, prevalentemente a. Basta salire al 
successivo orbitale libero per ritrovare il carattere tt: dunque, il 40(LUMO+3) rimpiazza il 39. 

Lo spazio attivo è allora dato dagli orbitali: 32,35,36,37,38,40. 

2.3 Metodi di scansione 

L'esigenza di sviluppare un metodo di scansione ad hoc per condurre il calcolo CASSCF delle 
superfici So ed Si sulle nostre griglie di punti nasce da due problemi. 

Il primo è che, come visto in sec.2.1.4, non è possibile definire un punto G Crea con le sole 
variabili della OJtairice Z , per via della dipendenza non lineare degli angoli esterni da P2 e P3 
e della natura collettiva di Z. Pur essendo implementati in Qaussian alcuni metodi di scansione 
(scan) capaci di generare griglie di punti, essi agiscono appunto sulle variabili della 9}tairice Z in 
modo indipendente e lineare dunque non soddisfano i nostri scopi. 
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2 Campionamento delle superfìci di energia potenziale 






Figura 2.11: Orbitali molecolari RHF disegnati con VMD [19] ad un valore di isosuperfìcie di 0.04. Il 36 è HOMO, il 37 
LUMO, seguendo la numerazione assegnata dal calcolo. 




Il secondo motivo è strettamente legato alle caratteristiche 
del metodo CASSCF, il quale necessita di preservare lo spazio 
attivo da un punto della griglia all'altro. Definito lo spazio at- 
tivo ad una geometria di riferimeno (vedi sezione precedente), 
la sua conservazione è ottenuta riutilizzando il set degli orbi- 
tali attivi di un dato punto come stima iniziale (guess) per un 
punto della griglia ad esso adiacente. 

Questo metodo è destinato a fallire qualora si cerchi di uti- 
lizzare gli stessi orbitali per due punti con geometrie molto 
diverse: Qaussian non riesce ad adattarli ed il calcolo diverge 
o restituisce risultati non fisici. Durante una scansione tradi- 
zionale, Qaussian legge gli orbitali del punto precedente, e passa al punto successivo modificando 
una variabile: se la scansione è monodimensionale basta scegliere step sufficientemente piccoli tra 
una geometria e l'altra e non si incontrano difficoltà. 

In una scansione bidimensionale, Qaussian esplorerà tutti i 
valori della prima variabile, quindi la azzererà per incrementare 
la seconda e ripartire daccapo con la prima. Quasi certamente 



Figura 2.12: Scansione tradizionale di 
Qaussian : tra un salto 
e l'altro le geometrie cam- 
biano troppo per il calcolo 
CASSCF 
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Figura 2.13: Possibile estensione di un 
metodo di scansione da 2 a 
3 dimensioni 



2.3 Metodi di scansione 



dopo questo step di azzeramento il calcolo CASSCF fallirà, 
perché gli orbitali memorizzati nel punto precedente saranno 
inconciliabili con il nuovo punto. 

Questi due problemi hanno reso necessario elaborare un pat- 
tern di scansione senza salti di variabile. Ci poniamo il pro- 
blema inizialmente in due dimensioni, per poi estenderlo a 3 
ed N come mostrato in fig.2.13. 

Ad esempio, lungo un pattern monosequenza ("worms" e "spirale" in fig.2.14) tutto un pia- 
no è esplorato da un singolo percorso, senza biforcazioni, che unisce ogni punto ad un solo 
punto precedente e ad un solo punto successivo. Il difetto di questo metodo sta nella scarsa 
parallelizzabilità. 
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r igura 2.14! Diversi tipi di scansione: monosequenza, a libro, ad onda 

Il pattern impiegato inizialmente per i nostri calcoli è quello a "libro": dal punto di partenza 
della scansione si diramano due assi di scansione lungo una variabile (la y in figura); di ogni 
punto appartenente a questi assi vengono memorizzati gli orbtiali al termine del calcolo. Comple- 
tato il primo stadio^ per ogni punto viene lanciata una scansione monodimensionale lungo x con 
l'algoritmo implementato in Qaussian . 

In tre dimensioni occorre conservare gli orbitali di tutti i punti di un intero piano, ed un terzo 
stadio esegue la scansione lungo la terza coordinata punto per punto. La scansione a libro ha 
una buona parallelizzabilità poiché il secondo stadio può cominciare mentre il primo deve ancora 
completarsi. 

Per raffinare e velocizzare ulteriormente il calcolo è successivamente sorta la necessità di con- 
servare gli orbitali di tutti i punti dello spazio. Con la scansione a libro questo non è possibile, 
perché l'ultimo stadio impiega la scansione monodimensionale di Qaussian che adatta gli stessi 
orbitali ad ogni punto successivo, perciò rimangono in memoria solo quelli dell'ultima geometria 
e non quelli di tutti i punti della griglia. 

Si è allora optato per una scansione "ad onda^", completamente gestita dal programma ester- 
no ReQuieM (app.A). La scansione si propaga senza una regola geometrica prestabilita: nel 
momento in cui è disponibile un processore per svolgere il calcolo, il gestore esplora lo spazio alla 
ricerca di un punto da completare che sia attiguo ad uno completato. Due punti sono attigui se 
hanno le stesse coordinate eccetto una, ed in quella differiscono di un solo step sulla griglia. 



^Così chiamata poiché il calcolo si propaga a sfera dall'origine del calcolo, e viene rallentato e deflesso dalle regioni 
dove richiede più tempo 
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2 Campionamento delle superfìci di energia potenziale 

In questo modo si soddisfa la condizione di minima variabilità tra geometria del punto da 
calcolare ed orbitali di partenza; si raggiunge la massima parallelizzabilità in quanto ogni punto 
completato rende possibile processare i 6 punti attigui; infine al termine si dispone di tutti i dati 
per calcoli futuri. 

L'ottima versatilità del linguaggio di programmazione con cui ReQuieM è scritto, Python 
[39], ha facilitato l'estensione del codice in modo da distribuire il calcolo in rete, verso tutte le 
workstation Linux (nodi) cui ci è stato garantito l'accesso presso l'ISTiTUTO per i Processi 
Chimico Fisici del CNR di Pisa, senza necessità di installare o configurare alcuna applicazione 
aggiuntiva sui nodi di lavoro (fatta eccezione per il preinstallato Qaussian ). 

2.4 Prima griglia: torsione e piramidalizzazioni 

La prima scansione è stata realizzata lungo la coordinata di torsione X e le due di piramidalizza- 
zione P2,3, sulla griglia di 1089 punti (XPP) data da tutte le combinazioni di: 

• X G {70,74,78,82,86,90,96,102,108} 

• P2^ {75,78,81,84,87,90,93,96,99,102,105} 

• Ps ^ {75,78,81,84,87,90,93,96,99,102,105} 

• Z=l (ferma alla CI) 

La griglia è centrata sull'intersezione conica poiché in questa regione la molecola ha maggior 
probabilità di decadere dallo stato eccitato a quello fondamentale: le caratteristiche della superficie 
di energia potenziale su cui i nuclei si muovono attorno alla CI è di importanza primaria per la 
dinamica di decadimento. 

Le energie calcolate da Qaussian in Hartree sono trasformate in kcal/mol e traslate dell'energia 
del minimo Trans: E = (E'^awssian — 402.08040^/i)*627.509A:ca//(7no/£^/j,); i dati sui due stati sono 
stati poi interpolati con Mathematica [62] in modo da ottenere una superficie continua. Il risultato 
della scansione sono due ipersuperfici di energia potenziale tridimensionali, non rappresentabili 
in grafico. In figura 2.15 sono riportate le PES in uno spazio cartesiano tridimensionale come 
funzione di due variabih per volta : (X,P3;P2=90°); (X,P2;P3=90°); (P3,P2; X=90°). 




Figura 2.15: Superfìci di energia potenziale della griglia XPP, in arancione per SI (superiore) e azzurro per SO (inferiore). 
Da sinistra in funzione di X e P3; P2=90° // (X e P2; P3=90°) // (P3 e P2; X=90°) 
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2.4 Prima griglia: torsione e piramidalizzazioni 

Con l'analisi delle PES così calcolate e con l'individuazione dei loro principali punti topologici 
lungo le tre coordinate scelte, cercheremo di rispondere alle due domande seguenti: 

1. Qual'è il ruolo di P2 e P3 sull'intersezione conica? È possibile trascurare una delle due, 
o almeno trattarla in modo più approssimato? Se l'effetto di una delle due piramidalizza- 
zioni fosse trascurabile il modello potrebbe essere ridotto da 4 a 3 dimensioni con notevole 
vantaggio computazionale ed interpretativo. 

2. Nel caso P2 e P3 siano entrambe rilevanti, sarebbe possibile introdurre una combinazione 
lineare di P2 e P3 che descriva adeguatamente la dinamica del sistema? Anche in questo caso 
avremmo una riduzione di dimensionalità con conseguente semplificazione del problema. 

2.4.1 Topologia dello stato eccitato 

Lo stato eccitato presenta un minimo di energia potenziale al punto X=90°,P2=90°,P3=90°, ov- 
vero il punto della griglia che si avvicina di più all'intersezione conica (89.53610°, 90.989663°, 
89.871629°). L'energia di Si (Es^) non diminuisce uniformemente attorno alla CI, bensì ha la for- 
ma di una valle con pareti di energia potenziale che salgono rapidamente eccetto in una direzione, 
lungo la quale l'energia cresce molto più lentamente formando un canale. 

Il canale è visibile in fig.2.15 in entrambi i grafici X-P, poiché corre lungo entrambe le coordinate. 
In figura 2.16 il canale è rappresentato dalle isosuperfici in cui l'energia eguaglia tre soglie diverse. 
Le zone a bassa energia si visualizzano più convenientemente disegnando dei grafici a curve di 




Figura 2.16: isosuperfici nelle tre coordinate in cui la E^^ è uguale a, dall'esterno al centro del volume, 
kcal/mol (la CI ha energia 58kcal/mol) 



/ 63 / 60 
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2 Campionamento delle superfìci di energia potenziale 

livello, dai quali si può anche stimare la pendenza della superficie: la figura 2.17 riporta curve di 
livello nello spazio (P3,P2) (cui ci riferiremo come PP) mantenendo X costante. 



EKcited state. CI geo 
at torsioni 00 



min: X=100, P2=92.2744. P3=8S.7608 



EKcited state. CI geo 
attorsioii89.5 





min: X=89.S. P2=91.02S2. P3=89.7201 



EKcited state. CI gì 
at torsion 80 




min: X=80, P2=89. 22G. P3=93.4076 



Figura 2.17: Curve di livello per E^^ in funzione di P3 (ascissa) e P2 (ordinate) a diversi valori di torsione (nel titolo) 



Dalle curve di livello si individua nettamente un minimo energetico unico nello spazio PP per 
ogni valore di X scelto. Per valori di torsione superiori a 90° questo minimo si trova nel quadrante 
(P3<90°,P2>90°); per X<90° invece il minimo si sposta nel quadrante simmetrico rispetto a 
(90,90), ovvero (P3>90°,P2<90°). 

I punti di minimo, calcolati per valori di torsione fissati tra 71° e 108° applicando un algoritmo 
di minimizzazione alla E^*^ interpolata, sono riportati in figura 2.18. 



Piramidalizations 
94 



m Excited State: Red=P3, Blue=P2 





Figura 2.18: A sinistra: Punti di minimo per E^^ calcolati nello spazio (P3,P2) a diversi valori di X. La torsione è in 
ascissa; la curva rossa indica l'angolo P3, quella blu il P2. A destra: rappresentazione degli stessi dati in 
forma di MEP; la sfumatura di colore indica la torsione, da sinistra X=71° a destra X=108° 

Tali punti non sono altro che il cammino di energia potenziale minima (MEP^) per l'isome- 
rizzazione sullo stato eccitato. Il MEP è una prima stima del percorso seguito dal CPP dopo 
l'assorbimento di un fotone, postulando che le energie termica e cinetica della molecola siano 
trascurabili rispetto a quella dell'eccitazione, nel corso di tutto il decadimento. Queste condizioni 



Dall'inglese Minumum Energy Path 
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2.4 Prima grìglia: torsione e piramidalizzazioni 

non sono realizzabili naturalmente, perché il processo è molto veloce quindi le energie cinetiche 
in gioco sono alte e la molecola può deviare ampiamente dal MEP. 

Nel tratto di torsione indagato la P3 del minimo su E^'^ passa da 93.5° a 86.5°, la P2 da 89° 
a 92°: la variabilità della prima piramidalizzazione è circa tripla rispetto a quella della seconda. 

Possiamo concludere che durante un movimento di isomerizzazione su Si la coordinata P2 varia 
molto meno rispetto alla P3 lungo il MEP. 

2.4.2 Topologia dello stato fondamentale 

Dalle curve di livello per lo stato fondamentale (fig.2.19) si osserva che l'energia di So (E^'q), 
similmente a quanto visto per E^'^, non diminuisce allo stesso modo in tutte le direzioni dando 
luogo ad una cresta nello spazio (X,P2,3) (vedi in fig.2.15) e ad uno o due minimi nello spazio PP, 

(P3,P2). 

Per 90° < x< 89° si notano due minimi nello spazio PP, come suggerito anche dalla leggera 
gobba visibile sulla isosuperficie a 57.5kcal/mol in fig.2.19. Non abbiamo dati abbastanza densi 
in quell'intervallo per stabilire se si tratti di un artefatto dovuto alla interpolazione: la griglia 
passa infatti solo per i punti X=86°, 90°, 96°. Nell'arco di quei tre step si registra però una forte 
discontinuità nelle coordinate di piramidalizzazione del minimo. 

Come già visto per Si, ricerchiamo il MEP minimizzando E^'q in funzione di P3 e P2 con X 
vincolata a valori tra 71° e 108°. L'algoritmo trova il minimo assoluto entro lo spazio richiesto 
perciò si ottiene un MEP globale, senza la biforcazione suggerita dalle fig.2.19. 

Il cammino di energia minima sullo stato fondamentale, analogamente a quanto già detto per 
lo stato eccitato, è una prima approssimazione del percorso che il CPP potrebbe intraprendere 
una volta decaduto dallo stato eccitato a quello fondamentale. 

Anche sullo stato fondamentale notiamo che il MEP percorre sulla coordinata di piramidaliz- 
zazione del pirrolinio molti più gradi rispetto a quella del pentene (tab.2.4). 

X=71° X=107° A 
P3 86.724° 93.380° 6.656° 
P2 91.012° 89.853° 1.159° 

Tabella 2.4: Spazio percorso dal MEP lungo le due coordinate di piramidalizzazione su E^q 

La variabilità del MEP rispetto a P3 è quasi sei volte superiore a quella di P2. 

2.4.3 Linee di quasi-degenerazione 

L'intersezione conica non è l'unico punto in cui è possibile il salto dallo stato eccitato a quello 
fondamentale: la probabilità di tale evento dipende infatti in una certa misura dalla differenza in 
energia potenziale tra gli stati. Introduciamo allora la superficie AE=FjSi-^So'- i minimi e le valli 
presenti su questa nuova funzione sono dette linee di quasi-degenerazione (LQD) tra i due stati. 
Le curve di livello (fig.2.20(a)) e le isosuperfici (fig.2.20(b)) indicano che esistono due linee di 
quasi-degenerazione dove la differenza AE è più bassa che nel resto dello spazio. Le isosuperfici 
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2 Campionamento delle superfìci di energia potenziale 



Ground state. CI geo 



Ground state. CI 



Ground state. CI geo 



Ground state. CI geo 




X=80, P2=91.3413, P3=86.411S 



P3 90 




Piramidalizations 
94 



Minima on Fundamental State: Red=pn. Blue=pc 



Cammino di energia minima su SO per 1': 





Figura 2.19: Sopra: Curve di livello per E^^ in funzione di P3 (ascissa) e P2 (ordinate) a diversi valori di torsione (nel 
titolo). Mezzo: Isosuperfici con E^^ uguale a (dall'esterno) 61.5, 57.5 e 53.5 kcal/mol. La zona di doppio 
minimo intorno a 89° è rintracciabile nella protuberanza sulla seconda isosuperfìcie. Sinistra: Punti di 
minimo per E^^ calcolati nello spazio (P3,P2) a diversi valori di X. Destra: MEP su E^^, la sfumatura di 
colore indica la torsione, da destra X=71° a sinistra X=108°. 
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2.4 Prima grìglia: torsione e piramidalizzazioni 

illustrano una struttura ad "x" o a "cromosoma" in cui le valli di A^ si incrociano nel punto di 
CI, intorno alla quale c'è la regione più profonda. 

I valori (P3,P2) di questi due pozzi (in fig.2.20(c)) sono stati calcolati mediante due minimizza- 
zioni nei piani PP ad X costante: la prima imponendo P2>90°, la seconda con P2<90°. In tal 
modo si sono ottenuti i due minimi locali di A^ per ogni X, che si collocano appunto in entrambe 
le possibili deviazioni dalla planarità dell'anello pentenico (P2 ^ 90°). 

La scarsa risoluzione dei dati non ci permette di dire se ad X=90° permangano due minimi 
separati o se invece i due canali coalescano. Trattandosi di una differenza di energie, quella che 
separa i pozzi non è tuttavia da considerarsi alla stregua di una barriera energetica con effetti 
dinamici. La molecola seguirà la dinamica determinata dall'energia potenziale dello stato su cui 
si trova. 

Per valutare dunque come CPP percorra queste linee di quasi-degenerazione le confrontiamo 
con i cammini di energia minima precedentemente calcolati per E^'^ ed E^'^ . I due MEP (al centro 
di fig.2.20(d)) sono molto più localizzati attorno a P3,2=90° rispetto alle LQD, e seguono un 
andamento trasversale in funzione di torsione e di P2 : si trova infatti a P2>90° per X<90°; a 
P2<90° per X>90°. La regione coperta è mediamente più vicina alla linea di quasi degenerazione 
P2>90°. 

Dalla figura 2.20(b) si nota come questa linea trasversale sia anche la più profonda, dove cioè 
la degenerazione è massima: le regioni AE< 2 kcal/mol compaiono infatti lungo il canale che va 
da (X=94°, P3=80°, P2=98°) a (X=83°, P3=105°, P2=81°) e non lungo l'altro. 

Per torsioni lontane da 90° le linee di quasi-degenerazione escono dalla nostra griglia per quanto 
riguarda P3,2' in particolare P3 esce per 94°<X<84°, mentre P2 per X<84° e solo nella linea 
P2<90°. A questo sono dovuti i molti punti accumulati all'inizio ed alla fine dei grafici delle linee. 

La valutazione riguardo la rilevanza delle coordinate di piramidalizzazione per descrivere le linee 
di quasi-degenerazione è dunque complessa. L'elevato gradiente di P3 almeno fin dove è definita 
la griglia, evidenziato in fig.2.21 a destra, lascia supporre che le linee di intersezione dipendano 
molto più strettamente da P3 che non da P2. Quest'ultimo svolge invece un ruolo più incisivo 
nella sfera di raggio 1-2° attorno dall'intersezione conica: l'isosuperfice più interna in 2.20(b) è 
infatti modellata principalmente lungo P2. Le altre isosuperfici invece coprono 20° su P2 e ben 
35°, ovvero tutti quelli disponibili in griglia, su P3. 

2.4.4 Regioni di rilassamento su So 

L'esito Cis o Trans del decadimento non radiativo da stato eccitato a fondamentale dipende dal 
momento posseduto dalla molecola quando si ritrova sullo stato fondamentale, e dal campo di 
forze in cui si deve muovere. Trascurando l'energia cinetica, una stima dell'esito del decadimento 
in un certo punto D=(X,P3,P2) la si ottiene combinando il segno della derivata prima rispetto ad 
X del potenziale E^'^ in D con quello di sinX: 

T = sign -^^\d sign sm X 

T > 0^ Cis ; T<0'w Trans 
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2 Campionamento delle superfìci di energia potenziale 



Diff between states. CI g&o 
at torsion 95 



Diff between states. CI geo 
at torsion 89.5 





Diff between states. CI geo 
ittiirsionSy 



. 75 




75 SO 85 90 95 100 105 75 80 85 90 95 100 105 75 SO 85 90 95 100 105 

mini: X=9S. P2=83.S9GS. P3=7S. mini: X=S9.S. P2=87.4S92, P3=87. mini: X=37. P2=87.. P3=9S.99S4 

min2: tor=95. p. pent=9G. 4139. p. pir=78. S985 min2 : tor=89.S, p. pent=92. 3647, p.pir=90. min2 : tor=87, p.pent=9G.. p.pir=99. 

(a) Curve di livello per AE in funzione di P3 (ascissa) e P2 (ordinate) a diversi valori di torsione (nel titolo). 




P3 



(b) Sinistra: Sovrapposizione delle isosuperfici per E^q edE^.^ . Destra: Isosuperfici per A^ uguale a (da quella esterna) 4, 3 
e 2 kcal/mol. 



Linee di intersezione 




(c) linee di quasi-degenerazione; la torsione aumenta da (d) Paragone MEP per So (marcato con O) ed Si 
sinistra a destra. Al centro, MEP sui due stati. (2.20(d)) (marcato con I) 

Figura 2.20: Analisi di AE 



40 



2.4 Prima griglia: torsione e piramidalizzazioni 




85 90 95 100 105 




80 85 90 95 100 105 



Figura 2.21: A sinistra: Andamento delle P2 delle due linee di intersezione, in funzione della torsione. A destra: per P3 

Se T > 0, allora il gradiente del potenziale è diretto verso l'isomero Cis (X=0°); altrimenti è 
diretto verso quello Trans. Il segno del seno occorre per localizzare quale minimo precede e quale 
succede a D lungo la coordinata periodica X. 

Calcolando T per ogni punto della griglia delineiamo due regioni di rilassamento RR in cui, 
applicando le stesse approssimazioni menzionate riguardo i MEP, la molecola decadrà con più 
probabilità verso un isomero o l'altro a seconda dell'T associato a quella regione. 

Intuitivamente ci aspetteremmo che dove il punto D ha torsione maggiore di 90°, la molecola 
decada sempre alla geometria Trans, e alla Cis per X<90°. Dal grafico tridimensionale delle RR in 
fig.2.22(a) si nota invece che non si incontrano lungo un piano perpendicolare all'asse X e passante 
per 90°. Lo spartiacque tra i due segni di T è una superficie quasi planare, perpendicolare al piano 
X,P3 ma inclinata rispetto al piano X,P2. 

In fig.2.22(b) sono calcolate le stesse regioni nello spazio PP applicando una minimizzazione 
locale a partire da ciascun punto D con X vincolata. Il calcolo segue il gradiente fino a giungere 
a Trans o Cis (T invece da informazione sul gradiente solo nel punto D). Il piano è colorato in 
funzione del punto di arrivo del minimizzatore, se Cis (rosso) o Trans (blu), e il risultato è del 
tutto simile a quanto ottenuto col metodo T. 

La dipendenza delle regioni di rilassamento dalla coordinata P2 è trascurabile poiché l'interfaccia 
tra regione Cis e Trans è approssimativamente parallela all'asse P2. Lo si vede in tre dimensioni 
(fig.2.22(a)) ed in due (fig.2.22(b)): ponendoci ad esempio ad X=92°, una piramidalizzazione di 
soli P3 = 84° modifica la RR (da Cis) in Trans, mentre lo stesso effetto non è ottenibile con la sola 
P2. Questo è possibile perché l'interfaccia tra le regioni sviluppa invece un certo angolo rispetto 
all'asse P3. 

In generale abbiamo che, per qualsiasi valore di torsione: 

• P3>90°, favorisce il rilassamento a Trans 

• P3<90°, quello a Cis. 

Più la torsione si allontana da 90°, più la piramidalizzazione deve essere consistente per invertire 
il rilassamento normalmente prevedibile da X (cioè, X>90°, trans; X<90°, cis). L'effetto di P3 
rimane valido per 85°<X<95°: dopodiché innanzitutto non sono disponibili dati, ma soprattutto 
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2 Campionamento delle superfìci di energia potenziale 




(a) 



Ground State Decay 
X=90. 



Ground State Decay 
tor=91,92 




7S Bi) m:^ ;•);) ^i5 lijo 

ELue-trar.5; Red- Ci s; V/hi te:- Locai 



105 75 80 85 90 95 100 

Elue=tran3; Red=CÌ3; White=Local 



105 75 so 85 90 95 100 105 

BLue=tran3; Red=CÌ3; White=LocaL 



(b) 

Figura 2.22: Regioni di rilassamento: in blu T < rilassa a Trans; in rosso T > rilassa a Cis 

è presumibile che il valore di P3 necessario per l'inversione esca dal campo di applicabilità della 
coordinata, a causa delle approssimazioni fatte in sez.2.1.7. 

2.4.5 Riduzione della dimensionalità 

Nelle sezioni precedenti abbiamo visto le seguenti proprietà delle PES su cui si muovono i nuclei 
del CFP: 
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2.5 Interazione torsione-piramidalizzazione 

(I) Il MEP sullo stato eccitato percorre lungo P3 una distanza tre volte maggiore rispetto a 
quella percorsa lungo P2. La dinamica sullo stato eccitato si muoverà dunque molto più su 
P3 che su P2 (sez. 2.4.1). 

(II) Sullo stato fondamentale l'importanza di P3 è ancora più marcata rispetto a P2, in quanto 
il MEP si sviluppa lungo la prima coordinata sei volte più di quanto si estenda lungo P2. 
La dinamica su So coinvolgerà P2 ancor meno di quella sullo stato eccitato (sez. 2. 4. 2). 

(Ili) Le linee di quasi-degenerazione percorrono entrambe le coordinate per tutta l'estensione 
concessa dalla griglia, incrociandosi ad "x" lungo P2. Alla CI la linea si sviluppa princi- 
plamente lungo P2. Tuttavia, la coordinata P3 delle linee varia in funzione di X molto più 
rapidamente di quella P2. Inoltre sovrapponendo alle linee le MEP sugli E^'^ ^ si ricava la 
linea dinamicamente preferita; questa dipende più fortemente da P3 che da P2 (sez. 2. 4. 3). 

(IV) La coordinata P2 è quasi ininfluente sulle regioni di rilassamento Cis/Trans, che dipendono 
invece fortemente da P3 (sez.2.5.1). 

In base a questi elementi, in particolare il IV, abbiamo deciso di trascurare la coordinata di 
piramidalizzazione del pentene P2 nella costruzione delle griglie future, onde diminuire la com- 
plessità ed il costo computazionale del calcolo. Riguardo al punto III, si accetta così di esplorare 
un piano perpendicolare a P2 tra le due linee di intersezione. 

2.5 Interazione torsione-piramidalizzazione 

Le regioni di rilassamento definite in sez.2.5.1 meritano ulteriore approfondimento, in quanto 
evidenziano in modo assai chiaro un'interazione tra la coordinata di torsione X e quella di pira- 
midalizzazione P3 capace di dirigere il rilassamento del CPP. La topologia di Esq indica come 
questa interazione torsione-piramidalizzazione (ITP) possa condurre la molecola fino all'isomero 
Cis oppure Trans in base al punto in cui essa decade Si ^ So, influenzando così la resa quantica 
della reazione. 

La pendenza verso Cis o Trans della superfice Esq per valori di torsione prossimi ai 90° di- 
pende notevolmente da P3: anche dove la torsione tra i due anelli fosse tale da favorire intuiti- 
vamente un rilassamento in una direzione, poniamo Trans, una piramidalizzazione del pirrolinio 
sufficientemente distorta potrebbe invece condurre la molecola nella direzione opposta, Cis. 

Nella presente sezione si interpreta questa interazione in base alla morfologia degli orbitali del 
sistema tt costituenti lo spazio attivo, arrivando ad una conclusione del tutto generale per questo 
tipo di legami coniugati e che correla, al momento solo qualitativamente, la resa quantica di 
reazione con gli orbitali molecolari. 

2.5.1 Morfologia degli orbitali molecolari 

Abbiamo cercato un'interpretazione chimica all'interazione torsione-piramidalizzazione (ITP) di- 
segnando i tre orbitali elettronici occupati dello spazio attivo in 9 conformazioni, su una griglia 
tale da esplorare le combinazioni per 3 valori di X=90 e 3 valori di P3 =90: 
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2 Campionamento delle superfìci di energia potenziale 

• Xg ( 88,90,92 ) 

• P3 e ( 80,90,100 ) 

Le densità elettroniche sono state calcolate da Qaussian mediante 9 calcoli CASSCF corrispon- 
denti; i grafici degli orbitali sono realizzati con VMD [19] a partire dai dati Qaussian . 

Dal calcolo CASSCF risultano 3 orbitali con differente occupazione elettronica e morfologi- 
camente molto differenti. In figura 2.23 sono riportati gli orbitali per l'intersezione conica, che 
condivide con le altre geometrie molte caratteristiche. L'orbitale maggiormente occupato (I) è 
sempre localizzato sull'anello pirrolinico, mentre gli altri due (II, III) hanno occupazioni simili e 
attraversano con propaggini il legame tra i due anelli, intrecciandosi ortogonalmente. 

Durante il rilassamento la molecola cercherà di riformare il legame tra i sistemi tt dei due anelli, 
che in CI è ridotto alle sole propaggini su 3,2 mostrate nelle figure 2.23(c),2.23(b). Quale direzione 
di torsione garantirà la massima sovrapposizione tra gli orbitali, ed il ripristino più rapdio del 
sistema tt planare? 

Se osserviamo l'orbitale II, 2.23(c) la direzione Trans produce un aumento di sovrapposizione, 
mentre la direzione Cis assottiglia ulteriormente le propaggini che collegano i sistemi tt dei due 
anelli. L'orbitale III, essendo ortogonale, indica invece la direzione opposta: è evidente che 
qualsiasi direzione prenda la torsione uno dei due orbitali debba divenire non legante a favore 
dell'altro. 

Gli orbitali II e III possono essere efficacemente visualizzati come bande elastiche che avvolgono 
l'asse torsionale e si ancorano a lati opposti dei due anelli; ciascun "elastico" impartisce così un 
momento torcente r lungo l'asse 3=2, favorendo una direzione di rilassamento. 

Le due bande elastiche, essendo ortogonali, producono coppie di forze opposte (coppia Cis Tc / 
coppia Trans Tt). L'ipotesi alla base della nostra interpretazione è che la coppia Cis o 
Trans impartita dagli orbitali II, III che legano i due sistemi tt sia in qualche modo 
proporzionale alla loro occupazione elettronica. Detto in termini più chimici, l'orbitale 
destinato a divenire non-legante durante la torsione è quello meno occupato; o ancora, la torsione 
favorisce il carattere legante dell'orbitale più occupato. 

L'intersezione conica ha una geometria ad alta simmetria dunque le occupazioni di II, III sono 
molto simili e così, secondo la nostra ipotesi, anche le coppie TcT^. Nelle pagine successive per- 
corriamo punto per punto la griglia osservando come variano le occupazioni e correlando le coppie 
generate da II,III^ con la regione di rilassamento in cui si trova ciascuna geometria. (Ricordando 
che P2=90°). 

La tabella 2.5 riassume, per ciascun punto della griglia: 

• La regione di rilassamento associata a quelle coordinate (ricordando che P3 è fisso a 90°), 
calcolata con il metodo T ). 

• Il momento torcente/coppia prodotto dall'orbitale più occupato dopo l'I che abbiamo tra- 
scurato: trans o cis 



"^L'orbitale I non è riportato poiché riguarda solo l'anello pirrolinico, senza essere interessato dalla torsione, ma è 
stato verificato in ogni punto. 
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2.5 Interazione torsione-piramidalizzazione 




(a) I, occupazione 1.99675, nessuna coppia 




(b) II, occupazione 1.92947, coppia Trans rt 




(e) II!, occupazione 1.92713, coppia Cis Tc 



r igura 2.23! I tre orbitali occupati nella geometria di intersezione conica. 
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2 Campionamento delle superfìci di energia potenziale 



Tabella riassuntiva degli orbitali ITP 


X 


P3 


Regione 


Coppia r II 


Coppia r III 


A occ. (II-III) 


Figura 


90° 


80° 


Cis 


Cis 


Trans 


0.00663 


2.24 




90° 
100° 


Trans 
Trans 


Trans 
Trans 


Cis 
Cis 


0.00234 
0.00995 


2.23 

2.25 


92° 


80° 


Cis 


Cis 


indet 


0.00610 


2.26 




90° 
100° 


Trans 
Trans 


Trans 
Trans 


Cis 
Cis 


0.00345 
0.01131 


2.27 
2.28 


88° 


80° 


Cis 


Cis 


Trans 


0.00708 


2.29 




90° 
100° 


Cis 
Trans 


indet 
Trans 


Cis 
Cis 


0.00205 
0.00902 


2.30 
2.31 



Tabella 2.5: Tabella riassuntiva degli orbitali II e III per l'interazione torsione-piramidalizzazione 

• La coppia generata dall'orbitale meno occupato: trans o cis. 

• La differenza tre le occupazioni II e III della matrice densità elettronica 

• il riferiento all'immmagine degli orbitali riprodotta nelle pagine successive 

Su tutti i punti della griglia l'ipotesi delle "bande elastiche" prevede correttamente la regione 
di rilassamento cui appartiene ogni geometria. L'unica incertezza si ha nel punto X=88°,P3=90° 
poiché l'orbitale II non presenta alcuna "banda elastica", e la direzione di rilassamento è invece 
espressa correttamente dall'orbitale III. 
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2.5 Interazione torsione-piramidalizzazione 




(a) II, occupazione 1.92745, coppia Cis Tc 



(b) III, occupazione 1.92082, coppia Trans rt 



Figura 2.24: x=90°, P3=80° 




(a) II, occupazione 1.92968 , coppia Trans rt 



(b) III, occupazione 1.91973, coppia Cis Tc 



Figura 2.25: x=90°, P3=ioo° 




(a) II, occupazione 1.92696, coppia Cis Tc 



(b) III, occupazione 1.92086, coppia indeterminata 



Figura 2.26: x=92°, P3=80° 
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2 Campionamento delle superfìci di energia potenziale 




(a) II, occupazione 1.92647, coppia Trans n (b) III, occupazione 1.923021, coppia Cis Tc 

Figura 2.27: x=92°, P3=90° 




(a) II, occupazione 1.93077, coppia Trans n 



(b) III, occupazione 1.91946, coppia Cis Tc 



Figura 2.28: x=92°, P3=ioo° 




(a) II, occupazione 1.92795 , coppia Cis Tc 



(b) III, occupazione 1.92087, coppia Trans n 



Figura 2.29: x=88°, P3=80° 
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2.5 Interazione torsione-piramidalizzazione 




(a) II, occupazione 1.92939, coppia indeterminata (b) III, occupazione 1.92734, coppia debolmente Cis Tc 

Figura 2.30: x=88°, P3=90° 




(a) II, occupazione 1.92910, coppia Trans rt 



(b) III, occupazione 1.92008, coppia Cis Tc 



Figura 2.31: x=88°, P3=ioo° 
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3 Diabatizzazione e forme funzionali 



L'analisi delle piramidalizzazioni svolta sulla griglia presentata nel capitolo precedente ha permesso 
di semplificare il modello quadri-dimensionale (4D) di CPP trascurando la piramidalizzazione del 
pentene P2, ritenuta meno significativa rispetto alla piramidalizzazione del pirrolinio P3, ottenendo 
così un modello in tre-dimensioni(vedi sez.2.4.5). 

In questo capitolo presentiamo il calcolo di una nuova griglia 3D (XPZ) ottenuta vincolando 
P2 e variando X,P3,Z; l'analisi delle PES ricavate permette l'individuazione delle caratteristiche 
topologiche principali. 

Si procede quindi alla diabatizzazione delle superfici ed alla ricerca di una forma funzionale per 
energia adiabatica e diabatica. 

3.1 Seconda griglia: Torsione, Piramidalizzazione, Stretching 

La seconda scansione è stata realizzata lungo la coordinata di torsione X, la piramidalizzazione P3 e 
lo stretching alternato Z, sulla griglia di 20x9x21=3780 punti (XPZ) data da tutte le combinazioni 
di: 

• X G {89.53601°; 10°i, < i < 18,i G Z}, ovvero da 0° a 180° con step di 10°, più il 
punto di CI esatta a 89.53601°. 

• P2 = 90.989663° (vincolata valore alla CI) 

• P3 ^ {89.871629° + 5°i, -4 < i < 4,i G Z}, da 69...° a 109...° con step di 5° 

• Z G { — , < i < 20,i G Z}, da a 2 con step di 0.1 

A parte le tre coordinate che formano la griglia, il resto della geometria è fissa a quella dell'inter- 
sezione conica. 

Poiché la griglia si estende fino ai punti cis e trans, abbiamo confrontato tutte le variabili in 
OJlaitice Z per le geometrie Trans e CI (in tab.3.1). Le differenze su variabili non incluse nelle 
nostre coordinate sono ampiamente trascurabili: le distanze di legame al di fuori del sistema tt 
variano un'ordine di grandezza in meno rispetto a quelle che prendono parte alla coordinata Z; le 
variazioni negli angoli e nei diedri, poi, sono tutte entro 1°. 

Per calcolare tutte le proprietà necessarie si è svolta una doppia scansione. Prima con Qaussian abbiamo 
determinato gli orbitali dello spazio attivo punto per punto e le distribuzioni di carica sui due 
anelli; quindi con-^LCAS 7, utilizzando gli orbitali calcolati da Qaussian , per calcolare i momenti 
dipolari su tutti e due gli stati elettronici con metodo RASSI. 
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3.1 Seconda griglia: Torsione, Piramidalizzazione, Stretching 



DISTANZE 


A CI- Trans 


Sistema n (Z) 




BIO 


0.03531 


B6 


-0.04188 


B2 


0.09031 


B5 


-0.05866 


B9 


0.02969 


Altri legami C-C 




B4 


-0.01778 


B7 


0.00358 


B8 


-0.00333 


Bll 


0.00041 


Legami C-H ibrid. sp2 


B16 


-0.00833 


B17 


-0.00617 


B18 


0.00611 


B19 


0.00491 


B12 


-0.00303 


B13 


-0.00269 


Altri legami C-H 




B14 


-0.00271 


B15 


-0.00235 


B20 


0.00447 


B21 


0.00434 


B22 


0.00693 


B23 


0.06077 



ANGOLI 


ACLTrans 


Pentene: P2 




A2 


0.98171 


A4 


0.98171 


A5 


0.98171 


Interni 




A8 


-0.49843 


A9 


-0.58804 


H legati al sìs.tt 




A16 


-0.14427 


A17 


1.35587 


Altri H 




A20 


-1.07539 


A21 


-1.00763 


A22 


-0.59537 


A23 


-0.22052 


Pirrolinio: P3 




A3 


-0.11856 


A7 


-0.11856 


A6 


-0.11856 


Interni 




AIO 
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DIEDRI 


ACLTrans° 


Torsione X 




D3 


-90.44776 


Esterni al Pentene 




D4 


-2.72148 


D5 


-0.64873 


Esterni al Pirrolinio 
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2.59320 


D7 


3.67011 


Altri 




D8 


-0.15804 


D9 


0.46230 


DIO 


0.03226 


Dll 


-0.04045 


D12 


-0.86258 


D13 


0.43899 


D14 


-0.01458 


D15 


-0.77337 


D16 


-0.13863 


D17 


0.5097 


D18 


0.15848 


D19 


-0.38077 


D20 


-1.13473 


D21 


1.25799 


D22 


-0.88802 


D23 


0.96130 



Tabella 3.1: Confronto tra tutte le variabili della VJlainct Z per le geometrie Trans e CI 

Entrambi i calcoli {Qaussian c-Wlcas ) hanno restituito delle PES: si è optato per utilizzare i 
risultati di-^LCAS in vista di una futura raffinazione del calcolo, con introduzione della correlazio- 
ne elettronica a livello CASPT2 (non ancora eseguita). Non c'è comunque differenza significativa 
tra le energie calcolate con-^LCAS e con Qaussian . 

3.1.1 Topologia delle PES 

La griglia XPZ condivide con la già studiata XPP la coordinata di torsione e quella di pirami- 
dalizzazione del pirrolinio. Ponendoci a Z=l ritroviamo perciò le stesse proprietà dipendenti da 
X,P3, principalmente i MEP, una linea di quasi-degenerazione tra gli stati e le due regioni di ri- 
lassamento sullo stato fondamentale. Come mostrato dalle curve di livello (fig.3.1(a),3.1(b)) e dai 
cammini di energia minima in funzione di Z, la topologia del piano X,P3 non viene essenzialmente 
modificata dal valore dello stetching. 

In figura 3.2 sono presentati i MEP sulle superfici di potenziale dei due stati elettronici adia- 
batici in funzione di X, ovvero per ogni valore di isomerizzazione abbiamo cercato il minimo sul 
piano P3,Z. Sullo stato fondamentale il minimo lungo la coordinata Z (fig.3.2(b)) si sposta, quasi 
linearmente, da Z=0 in Trans a Z=l in CI, per poi tornare a Z=0 in Cis. Viceversa su E^'^ il 
minimo si trova a Z=1.5 negli isomeri Trans e Cis, passando per Z=l nella CI. 

Il minimo lungo la coordinata di piramidalizzazione invece (fig.3.2(a)) si sposta da 90° di Trans, 
fino a P3= 95° pochi gradi di torsione prima della CI, saltando a P3=85° dopo la CI per poi risalire 
fino a P3=90° nella geometria Cis. 
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3 Diabatizzazione e forme funzionali 

Sullo stato eccitato Si, viceversa, il MEP decresce lungo la piramidalizzazione dai P3=90° della 
geometria Trans fino a P3=85° pre-Ci, quindi aumenta a P3=95° dopo la CI (più gradualmente 
rispetto al MEP su E^'^) per tornare a P3=90° in Cis. 
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(a) Curve di livello E^^ in funzione di X (ascissa) e P3 (ordinata), a cinque valori di Z 0, 0.5, 1, 1.5, 2 
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(b) Curve di livello E^^ in funzione di X e P3 
Figura 3.1: Sostanziale invarianza dell'ITP rispetto allo stretching Z 






(a) Visuale X,P3 dei MEP 



(b) Visuale X,Z dei MEP 



(e) Visuale P3,Z dei MEP 



r igura 3.2: Diverse visuali dei cammini di energia minima (MEP). Lo stato fondamentale è in sfumature di blu funzioni 
della torsione, dal chiaro (X=180°, trans) allo scuro (X=0°, cis). Quello eccitato va dal giallo (X=180°) al 
rosso (X=0°). 

Le superfici di energia potenziale nello spazio XPZ presentano solamente una linea di quasi- 
degenerazione (fig.3.3(a)), dipendente principalmente da X,P3 mentre la coordinata Z ha l'effetto 
modulare la differenza di energia (fig.3.3(b)). Il contributo di Z diventa fondamentale nel punto di 
intersezione conica: in fig.3.4 sono riportate in grafico le due superfici ad X,P3 nel punto di CI, e si 
nota come a Z=l le energie formino una cuspide in cui la loro differenza diminuisce rapidamente. 

3.2 Funzione d'onda e diabatizzazione 

Come sottolineato nella sezione precedente, nel punto di intersezione conica le superfici di energia 
potenziale calcolate sui due stati elettronici presentano una discontinuità, percui variano rapi- 
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3.2 Funzione d^onda e diabatizzazione 





(a) Curvedilivello A^(X,P3) aZ=l. Nella regione (b) Curve di livello AE{X,Z) a P3=89.9. Nella 
bianca A^=3.5 regione bianca A^=3.5 

Figura 3.3: Linee di quasi degenerazione nello spazio XPZ 
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Figura 3.4: Energie dei sue stati in funzione di Z (X=90°, P3=90°); in rosso E^o, nero E^i. In Z=l raggiungono con una 
cuspide la massima degenerazione. A destra: differenze tra le energie. 



damente in funzione delle coordinate nucleari XPZ. Questo pone difficoltà nel momento in cui 
si vogliano utilizzare le PES per condurre calcoli di dinamica, poiché nel punto di singolarità 
gli accoppiamenti non-adiabatici (derivativi) tra i due stati tendono ad infinito e rendono così 
impossibile il calcolo. 

Il concetto di stati elettronici diabatici ha un ruolo centrale nel trattamento dinamico di processi 
non-adiabatici [23]. Uno stato diabatico è definito in modo da non soffrire della forte dipenden- 
za dalle coordinate nucleari in prossimità di una intersezione conica: gli accoppiamenti cinetici 
non-adiabatici tra gli stati adiabatici, singolari attorno alla CI, vengono invece sostituiti da un 
accoppiamento diabatico di tipo potenziale tra i due stati diabatici, che varia in maniera smooth 
con le coordinate nucleari, anche nella regione della CI. 

Stati rigorosamente diabatici generalmente non possono essere definiti per molecole poliatomi- 
che come il CPP; ciò nonostante è possibile approssimare stati diabatici accettabili per lo studio 
della dinamica, detti stati quasi- diabatici. 
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3 Diabatizzazione e forme funzionali 

Gli stati quasi-diabatici possono essere calcolati con tre classi di strategie, a seconda dell'accu- 
ratezza e del costo computazionale: il metodo più elaborato prevede il calcolo degli accoppiamenti 
derivativi nelle coordinate rilevanti. 

Poiché gli accoppiamenti non-adiabatici diventano estremamente complessi da ricavare per gran- 
di sistemi molecolari, il secondo metodo sfrutta invece esclusivamente l'energia adiabatica ed i 
gradienti nella differenza tra le energie attrono alla CI. In tal modo non è necessario svolgere 
dispendiosi calcoli di quantistici aggiuntivi, ma si sfruttano esclusivamente le energie che già si 
devono già conoscere prima di qualsiasi calcolo dinamico. 

Una via di mezzo, in termini di affidabilità dei risultati e costo, è considerare informazioni più 
facilmente calcolabili riguardanti le funzioni d'onda elettroniche: proprietà in grado di indicare il 
mescolamento tra gli stati elettronici in funzione della geometria, e utilizzabili perciò per definire 
una trasformazione che possa regolarizzarli in stati quasi-diabatici, dipendenti in modo più smooth 
dalla geometria. 

Il CPP è una specie carica positivamente, e la distribuzione di questa carica sulla molecola 
dipende dalla funzione d'onda: perciò la carica potrebbe essere una buona proprietà per avere 
informazioni sullo stato elettronico del CPP. Per stabilirlo, abbiamo indagato le proprietà della 
carica nello spazio XPZ: per ogni punto della griglia abbiamo sommato le cariche di tutti gli atomi 
di ciascun anello in due numeri, ChPi (carica sul pirrolinio) e ChPe (carica sul pentene). Le figure 
3.5 illustrano l'andamento delle cariche in funzione delle tre coordinate. 

Allo stato fondamentale e nelle geometrie di bassa energia, la carica è prevalentemente sull'a- 
nello pirrolinico, ChPi>ChPe. In posizioni più prossime alla intersezione conica la differenza di 
carica tra i due anelli diminuisce gradualmente, e nelle immediate vicinanze alla CI ci sono forti 
discontinuità lungo tutte e tre le coordinate. Distorsioni lungo P3 e soprattutto Z(>1) porta- 
no la carica sull'anello pentenico, ChPe>ChPi, come mostrato in 3.5(a). L'andamento su Si è 
complementare al precedente: il passaggio di stato in una certa geometria, dunque, trasferisce la 
carica da un anello all'altro, eccetto che in quelle conformazioni dov'è distribuita metà e metà tra 
pentene e pirrolinio. 

La frazione di carica sui due anelli è una proprietà alquanto rudimentale su cui basarci per 
diabatizzare le PES: ad esempio, essendo una somma bruta non tiene conto dell'atomo sul quale 
si trova effettivamente la carica. 

Il momento dipolare è un operatore monoelettronico facilmente calcolabile per tutti i punti della 
nostra griglia, impiegando il metodo RASSI di-WLCAS 7. Dal calcolo si ottengono tre matrici i cui 
elementi rappresentano le componenti dei momenti dipolari sui due stati ed il momento dipolare 
di transizione tra i due, per ciascuna coordinata cartesiana. 



So ,,So ,,So\ 



I^So = (Mr,/^y%M. 



So So ,,So) 



i^s, = {fi^^°,i^;°,fit' 



tr ,,tr ,.tr\ 



f^tr — [f^x ^ l^y 5 /^, 
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3.2 Funzione d^onda e diabatizzazione 





70 80 90 100 no 70 80 90 100 

(a) Coordinata P3 in ascissa e Z in ordinata. A sinistra So, destra Si 
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(b) Coordinata X in ascissa e Z in ordinata. A sinistra So, destra Si 
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(e) Coordinata X in ascissa e P3 in ordinata. A sinistra So, mezzo Si. A destra: A^ nelle stesse coordinate; la linea di 
quasi degenerazione è sovrapponibile a quella lungo la quale le cariche si scambiano tra gli anelli e raggiungono la massima 
separazione (aree chiare) 

r Igura 3.0! Distribuzione di cariche sui due anelli. Le curve di livello indicano la differenza tra la carica totalizzata 
dall'anello pirrolinico e quello pentenico. In sfumature blu differenze positive, con più carica localizzata sul 
pirrolinio. In sfumature rosse differenze negative, con più carica spostata sul pentene. In colonna a destra le 
cariche su So, a sinistra su Si. Le sfumature più chiare indicano massima separazione. 
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3 Diabatizzazione e forme funzionali 




Figura 3.6: Distribuzione di carica tra i due anelli intorno alla CI, in funzione di una coordinata per volta. Curve blu e 
ciano: carica sul pirrolinio rispettivamente su So ed Si; curve rossa e arancione: carica sul pentene su So ed 

Si. 

•^LCAS infatti calcola queste proprietà in un sistema di riferimento cartesiano, quindi le tre 
coordinate x,y,z corrispondono agli assi cartesiani arbitrariamente scelti dal programma. Abbiamo 
dunque fissato tre assi nel nostro sistema di coordinate interne 3.7, con origine in C2: 

• Aa; è il versore definito dalla differenza vettoriale C3-C2 

• A^ è il versore definito da X1-C2 

• A^; è perpendicolare ai due precedenti, diretto nel semispazio contenente C5 




Figura 3.7: I tre assi cartesiani 

Per ogni punto della griglia si è calcolata la trasformazione necessaria per riportare C2 in 
origine, C3 lungo l'asse X e XI lungo l'asse Z. La matrice di rotazione torvata è invariabilmente 
la seguente: 



R = 



/ 9.99850824e-i -5.60541442e-^ -1.72631440e-2 
9.99473253e-i -3.245331086-^ 

\ 1.72722421e-2 3.24484696e-2 9.993241556-^ 



Applicando la rotazione R alle tre componenti fiSo , fJ^Si -, l^tr otteniamo dei momenti dipolari 
definiti rispetto a degli assi significativi nelle nostre coordinate interne. 

Le tre componeti di jisq-, l^Si^ l-^tr sono riportate nei grafici in fig.3.8 come funzioni delle tre 
coordinate di reazione XPZ. 

Concentriamoci sulle componenti M^^ della matrice momento di dipolo, ovvero quella diretta 
lungo l'asse A^^ (il primo grafico di ciascuna fig.3.8(a),3.8(b),3.8(c)). 



56 



3.2 Funzione d^onda e diabatizzazione 




(a) Le tre componenti di fiSo in funzione (da sinistra) di Z,P,X mantenendo le altre due coordinate fisse alla CI 




rO 1,5 2,0 -2- 

(b) Le tre componenti di fiSi in funzione (da sinistra) di Z,P,X mantenendo le altre due coordinate fìsse alla CI 




(e) Le tre componenti di fitr in funzione (da sinistra) di Z,P,X mantenendo le altre due coordinate fìsse alla CI. Il 
segno del momento di transizione è calcolato arbitrariamente in Molcas. 

Figura 3.8: Componenti dei momenti dipolari sui due stati in funzione di Z, P, Z. Colori: /2x in blu, /2y in viola, /2z in 
ocra. 
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3 Diabatizzazione e forme funzionali 

• L'asse A^^ indica la posizione della carica; il trasferimento da un anello all'altro percorre 
infatti proprio la direzione di quest'asse. 

• M^:, diretto lungo A^^, è la componente che varia più drasticamente, in funzione di tutte e tre 
le coordinate; questo è particolarmente evidente per la componente extradiagonale momento 
dipolare di transizione (trascurando i cambiamenti di segno arbitrari, vedi sez.3.3.1) 

• La rottura del legame 3,2 si ripercuote sul sistema tt ridistribuendo elettroni, carica e modi- 
ficando la funzione d'onda. Il sistema coniugato si sviluppa principalmente lungo l'asse A^^, 
e le variazioni di distanze nucleari possono ben essere approssimate a differenze di posizione 
lungo quell'asse. 

• In funzione di Z, M^^ riproduce qualitativamente l'andamento delle cariche: ciò conferma 
che Z sia la coordinata principale lungo la quale gli stati elettronici adiabatici si trasformano 
l'uno nell'altro e che M^^ ne indichi la maggior parte dell'effetto sulla funzione d'onda. 

Per questi motivi abbiamo scelto M^^ come proprietà della funzione d'onda indicatrice dello 
stato elettronico. 

La diabatizzazione consiste ora nel cercare una trasformazione che diagonalizzi M^^: i nuovi stati 
diabatici dovranno avere l'uno momento dipolare fortemente direzionalizzato verso il pirrolinio, 
l'altro verso il pentene, e momento dipolare di transizione nullo. Ragionando in termini di cariche, 
nella base diabatica il pirrolino ospiterà la carica positiva in uno stato ed il pentene nell'altro. 

Per ogni punto della griglia calcoliamo la matrice di trasformazione D che diagonalizza M^^, 
restituendo una matrice diagonale i cui elementi /i^, /i^ costituiscono i momenti dipolari diabatici. 



D^M^D ={ ^: (3.1) 




Il risultato della diabatizzazione è illustrato nelle figure 3.9: il metodo è efficace poiché i momenti 
dipolari diabatici sono ben separati nei due stati, ed in tutte le loro componenti - nonostante si 
sia diagonalizzata solo la componente M^^ del momento. 

Applicando ora la trasformazione D alla matrice Hamiltoniana otteniamo le energie corrispon- 
denti agli stati diabatici (Vii, V22) e l'accoppiamento tra i due stati (Vi2). 



D^HD = I ^'' ^'' ] (3.2) 

V V21 V22 J 

Le superfici diabatiche sono presentate in fig.3.10; si intersecano nelle zone in cui abbiamo 
visto scambiarsi le cariche ed i momenti dipolari. In funzione di Z (fig.3.10(b)) è particolarmente 
evidente la differenza tra stati adiabatici e diabatici: i primi presentano una cuspide all'intersezione 
conica per Z=l, senza attraversarsi nel rispetto della non-crossing rule. Le energie diabartiche si 
intersecano invece procedendo in modo smooth il loro andamento parabolico. 
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3.2 Funzione d^onda e diabatizzazione 
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(a) Le tre componenti (da sinistra X, Y, Z) di /isq in blu e /iSi in viola, in funzione di Z 



Componente X, diabatissata 



C omp onente Y, diab atis s ata 



Componente Z, diabatiEEata 



IrO . ^ 1.5 











0.1 












0.5 


1.0 1.5 2.0 


-0.1 


—^ 


"*-*-— ^ 


--^-^-. 


-0.2 









0.5 1.0 



1.5 2.0 



(b) Momento diabatizzato: le tre componenti (da sinistra X, Y, Z) di fia in blu e /if, in viola, in funzione di Z 

Componente Z 



Componente X 



Componente Y 
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(c) Le tre componenti (da sinistra X, Y, Z) di /iSo in blu e /isi in viola, in funzione di X 

Componente Z, diabatissata 



Componente X, diabatissata 



C omp onente Y, diab atis s ata 



■ 100 . . . 150 



.\ 




(d) Momento diabatizzato: le tre componenti (da sinistra X, Y, Z) di /la in blu e /i^ in viola, in funzione di X 
r Igura 3.9! Confronto tra gli elementi diagonali della matrice Mx con la proiezione sugli stati diabatici, ovvero fii e fi2 
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3 Diabatizzazione e forme funzionali 




(a) Superfici diabatiche in funzione di due coordinate per volta, la terza al valore della CI: da sinistra, (P3,Z), (X,P3),(X,Z) 




(b) L'energia adiabatica (sx) e diabatica (dx) in funzione di Z, (X e P3 fìsse alla CI) 
Figura 3.10: Le superfìci di energia potenziale negli stati diabatici 

Confronto con il modello di diabatizzazione di Stock 

Lo schema di diabatizzazione qui impiegato è sostanzialmente diverso da quanto scritto in let- 
teratura riguardo a molecole simili. In particolare Stock [18] ha proposto per il cromoforo della 
rodopsina un modello a due coordinate di reazione di cui una di torsione ed una di stretching, 
paragonabili alle nostre X e Z. 

Nella sua interpretazione la coordinata X è la quella reattiva responsabile della isomerazzazione 
e Z è una coordinata che accoppia i due stati elettronici. Nel costruire la matrice diabatica 
Stock pone un termine di accoppiamento linearmente dipendente da Z, ed i termini diagonali 
essenzialmente determinati da X. Questo da luogo a due superfici diabatiche che si intersecano per 
tutta la linea X=90° indipendentemente da Z; uno stato diabatico rappresenta così la geometria 
Trans e l'altro quella Cis. Nel nostro modello invece l'accoppiamento tra le energie diabatiche 
dipende quasi esclusivamente dalla torsione: la piramidalizzazione influisce 10 volte meno e lo 
stretching ha un effetto di due ordini di grandezza più piccolo ancora. Lo stretching, non la 
torsione, è la coordinata lungo la quale gli stati diabatici si scambiano, infatti del momento dipolare 
è stata scelta la componente M^ diretta lungo l'asse che maggiormente interessa lo stretching. 

Secondo la nomenclatura dei modelli vibronici di Cederbaum et cow. [63] la coordinata Z è quella 
di tuning: lungo Z le superfici adiabatiche si scambiano carattere elettronico ma non c'è quasi 
accoppiamento tra le diabatiche. Invece, quella torsionale è di accoppiamento: le due superfici 
adiabatiche non cambiano carattere elettronico ma c'è un forte accoppiamento diabatico. 

Riteniamo che il modello impiegato abbia maggior senso fisico rispetto a quello proposto 
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3.2 Funzione d^onda e diabatizzazione 




Figura 3.11: Diabatizzazione dei nostri stati adiabatici seguendo il modello proposto da Stock. In funzione di X,Z a 
P3=89.9° 



Debye 




0.5 1.0 



Figura 3.12: Accoppiamento in funzione delle tre coordinate (da sinistra: Z, P, X). A sinistra, l'accoppiamento è talmente 
piccolo da essere instabile numericamente. 
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3 Diabatizzazione e forme funzionali 

da Stock, perché descrive Cis e Trans con lo stesso stato diabatico perché, dal punto di vista 
elettronico, le due geometrie sono molto simili. 

3.3 Funzione analitica per gli stati adiabatici 

Dalla trasformazione eq.3.2 otteniamo due superfici di energia potenziale diabatica Vii,V22 (fig.3.10(a)) 
ed una superficie di accoppiamento V12. Al fine di sfruttare queste superfici in programmi di di- 
namica semiclassica e quantistica, nonché in campi di forza per la meccanica molecolare, occorre 
trovare una forma funzionale computazionalmente leggera che possa restituire tutte le proprietà 
necessarie a quei tipi di calcolo, quali: 

• Energie adiabatiche 

• Energie Diabatiche ed accoppiamento 

• Derivate parziali delle grandezze precedenti rispetto ciascuna delle tre coordinate Crea 

La via più diretta per ottenere questo tipo di forma funzionale è quella del fit lineare delle tre 
superfici, che possiamo calcolare punto per punto in tutta la griglia. Dai tre fit è facile ottenere 
le derivate diabatiche analiticamente; applicando la trasformazione inversa alla diabatizzazione 
possiamo calcolare le energie adiabatiche; infine dalla trasformazione inversa in forma derivata si 
ottengono le derivate delle energie adiabatiche. 

Come accennato in fig.3.8(c), il momento di transizione calcolato con -^lcas ha un segno 
arbitrario in ogni punto perché i calcoli sono svolti indipendentemente. Esso è definito in base al 
segno delle funzioni d'onda perciò, se questa informazione fosse preservata da un punto all'altro, 
non dovrebbe avere il comportamento random che invece osserviamo. 

Il segno di fitr non ha alcun effetto sulle energie diabatiche, ma determina invece il segno 
dell'accoppiamento (Acc). Se ne ottiene una superficie assai irregolare, impossibile da fittare 
correttamente (fig.3.13 a sinistra). 



Debye 




Figura 3.13: Accoppiamento in funzione di X calcolato a Z=0.8, P3=99.9. A sinistra e nel mezzo, col segno datogli da 
Molcas; a destra la stessa linea con il segno regolarizzato. 

Scegliamo il segno dell'accoppiamento di ciascun punto sulla griglia in modo da formare la su- 
perficie più regolare possibile, per restituire continuità al segno della funzione d'onda e semplificare 
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3.3 Funzione analitica per gli stati adiabatici 

al massimo il fit. Il metodo impiegato parte dalla considerazione che l'accoppiamento deve essere, 
al prim'ordine, una retta passante per zero nel punto di intersezione conica [32,49]. Percorriamo 
dunque l'Acc lungo una linea di torsione, avendo fermato le altre coordinate, e identifichiamo il 
punto di Acc minimo, Xmin. Imponiamo quel punto come ad Acc positivo assieme a tutti quelli 
aventi X>Xmin, e scegliamo il segno negativo per tutti gli Acc di punti con X<Xmin. Ne risulta 
un piano dal fit semplice ed efficace, in fig.3.14. 

La regolarizzazione sull'accoppiamento diabatico per imporre il comportamento atteso restitui- 
sce continuità al momento di transizione extradiagonale fi^r su base adiabatica, come dimostrato 
da Santoro. 





Figura 3.14: Accoppiamento regolarizzato: superficie planare di facile fit 



Il fit è stato svolto con il programma Mathematica [62] ricercando il minimo errore per una 
combinazione lineare, con coefficienti Xn,m,qi della funzione F: 

Wx = {cos(nX) - cos(n89. 536010)}, < n < 6 
Fp3 = {(P3 - 89.871629)^}, < n < 4 
F^ = {(Z-1)^}, 0<n<4 

6 4 4 

F = ^ ^ ^ Xn^rnAcos{nX) - cos(n89.536010))(P3 - 89.871629)^(Z - 1)H 

n=l m=0 q=0 
4 4 

Yl E Vm,,(i^3 - 89.871629)-(Z - 1)^ 

m=0 q=0 

Il fit conterrà dunque termini dipendenti da una o due coordinate contemporaneamente, per tener 
conto degli accoppiamenti tra le coordinate (ad esempio l'ITP in sez.2.5). 

Per il fit dell'accoppiamento si sono utilizzate solo le prime due funzioni di P3 e Z ed i primi 
4 coseni di X. Nei paragrafi successivi sono riportati i valori dei coefficienti per le funzioni sem- 
plificate: le potenze sono espanse e le potenze dei punti di equilibrio aggiunte al termine noto, 
riducendo F alla forma : 

6 4 4 
n=0 m=0 q=0 
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3 Diabatizzazione e forme funzionali 



Coefficienti per la superficie diabatica V 



11- 



Termine noto: k = 3963.32715674951 



Funzioni cos(nX): 


n 


coeff. 


1 


+217.115382 


2 


+ 1332.778533 


3 


-215.800246 


4 


-391.217098 


5 


+8.026211 


6 


-72.063152 



Funzioni (P^): 


m 


coeff. 


1 


164.220929 


2 


+2.630624 


3 


0.019107 


4 


+0.000053 



Funzioni (Z^): 
q coeff. 

1 -4440.787137 

2 915.756070 





Funzioni miste cos{nX){P^): 




n i, m - 


-^ 1 


2 


3 


4 


1 


-6.93972e+00 


6.78326e-02 


-1.36717e-04 


-6.65753e-07 


2 


-6.21455e+01 


1.06548e+00 


-8.05571e-03 


2.26661e-05 


3 


7.63456e+00 


-9.62699e-02 


4.90385e-04 


-7.57488e-07 


4 


1.73302e+01 


-2.89677e-01 


2.15487e-03 


-6.01461e-06 


5 


-1.43231e+00 


4.16536e-02 


-4.40109e-04 


1.58748e-06 


6 


3.71387e+00 


-7.00942e-02 


5.77522e-04 


-1.76025e-06 





Funzioni miste cos(nX)(Z"^): 


ni, 


q- 1 


2 


1 


4.77253e+00 


-2.72607e+00 


2 


3.40748e+01 


-2.44841e+01 


3 


3.20348e+00 


-1.59079e+00 


4 


5.76866e+00 


-3.91309e+00 


5 


6.55727e-03 


9.19716e-01 


6 


-1.01410e+00 


2.57722e+00 



Funzioni miste cos(P3"^)(Z"^): 
m i, q^ 1 2 



1.98270e+02 -3.31695e+01 

-3.33780e+00 5.38040e-01 

2.49256e-02 -3.98825e-03 

-6.96560e-05 1.13765e-05 



Coefficienti per la superficie diabatica V 



22- 



Termine noto: k = 3998.99093051 



Funzioni cos(nX): 


n 


coeff. 


1 


297.820058 


2 


175.1292169 


3 


-75.8710489 


4 


87.037889 


5 


173.921461 


6 


139.160875 



Funzioni (P^): 
1 coeff. 



-165.298438 

2.649808 

-0.019189 

5.3033413e-05 



Funzioni (Z^): 
q coeff. 

1 -4440.787137 

2 915.756070 





Funzioni 


miste cos(nX)(P3"^): 




n i, m 3 


-^ 1 


2 


3 


4 


1 


-1.14851e+01 


1.65009e-01 


-1.04123e-03 


2.42970e-06 


2 


-6.20786e+00 


9.14390e-02 


-6.56617e-04 


1.95312e-06 


3 


2.57821e+00 


-2.71020e-02 


6.91250e-05 


2.11920e-07 


4 


-3.91334e+00 


6.66198e-02 


-5.12016e-04 


1.49880e-06 


5 


-7.80964e+00 


1.31896e-01 


-9.92634e-04 


2.80838e-06 


6 


-6.19493e+00 


1.02776e-01 


-7.54977e-04 


2.07333e-06 
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Funzioni miste cos{nX)(Z^): 


n Xi, q ^ 


1 


2 


1 


-5.23070e-01 


-2.49707e+00 


2 


-1.43252e+01 


3.88926e+00 


3 


4.93237e-01 


-3.05146e+00 


4 


3.04809e+00 


-5.93955e+00 


5 


8.63719e-01 


-1.96025e+00 


6 


2.09983e+00 


-3.47894e+00 



3.3 Funzione analitica per gli stati adiabatici 



Funzioni miste cos(P3"^)(Z"^): 
m i, q ^ 1 2 



-2.78150e+01 9.90596e+01 

5.74679e-01 -1.77986e+00 

-4.76942e-03 1.37059e-02 

1.38016e-05 -3.83091e-05 



Coefficienti per la superficie di accoppiamento V 



Termine noto: k, ■ 



-16.7342746111 



Funzioni cos(nX): 
n coeff. 



-87.1362777844 
12.5694178958 
8.49199498388 
6.17103219602 



Funzioni (P^)'. 
m coeff. 



1 0.0974926541984 

2 0.00103044620733 



12- 



Funzioni (Z^): 
q coeff. 



-12.3455722879 
6.93514484091 



Funzioni miste cos{nX){ 


p^y. 


n Xi, m Ps ^ 1 


2 


1 7.43210e-01 


-4.09328e-03 


2 -1.10422e-01 


-3.02076e-04 


3 -9.75932e-02 


5.45871e-04 


4 -7.87783e-02 


1.57733e-04 



Funzioni miste cos(nX)(Z"^): 


n Xj, q Z^ 1 


2 


1 2.17239e+01 


2.52427e+01 


2 -3.69440e+00 


3.64576e+00 


3 -1.48621e+01 


1.44346e+01 


4 -2.11647e+00 


1.97404e+00 



Funzioni miste cos(P3"^)(Z"^): 
mPi, qZ^ 1 2 



2.71814e-01 3.98422e-02 
-1.75528e-03 -1.00075e-03 



Stima dell'errore 

L'errore associato alle tre superfici diabatiche, ora in forma analitica, si calcola prendendone i 
valori per ogni punto della griglia e ritrasfromando le energie in adiabatiche. Tra i valori adiabatici 
provenienti da-^LCAS e quelli calcolati dalle diabatiche fittate si prende il modulo della differenza 
come stima dell'errore. 

Sullo stato fondamentale l'errore medio ammonta a 0.287787 kcal/(mol/punto), su quello ec- 
citato a 0.398211 kcal/(mol/punto). L'andamento dell'errore è mostrato in fig.3.15. Si nota un 
errore maggiore intorno alla CI, che comunque rimane entro l-2kcal/mol. 



3.3.1 Implementazione in Fortran: E 



pisa 



Le forme funzionali presentate sono state implementate in un programma Fortran chiamato 
"Episa", listato in tab.3.3. Le routine ivi contenute sono state utilizzate nell'ambito della collabo- 
razione con Manuel Domer ed il prof. Daniele Passerone e per il calcolo di dinamica semiclassica 
affrontato nel prossimo capitolo. 

Il listato comprende solo il controllo da riga di comando e la subroutine episa in cui dalle 
energie diabatiche e dall'accoppiamento si ricavano le energie adiabatiche, con relative derivate. 
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3 Diabatizzazione e forme funzionali 




Figura 3.15: Distribuzione dell'errore dovuto al fit delle superfici diabatiche, calcolato rispetto a quelle adiabatiche. 



Nel resto del programma sono definite le funzioni F per ciascuna superficie diabatica dial, dia2 
e per l'accoppiamento coup; quindi sono definite le derivate di queste rispetto alle tre coordinate 
x,y,z: dialdx, dia2dx, coupdx, dialdy, dia2dy, coupdy, dialdz, dia2dz, coupdz. 

Il programma può essere compilato come eseguibile da riga di comando. In questa modalità 
Episa legge una riga dallo standard input contenente i tre valori X, P3, Z separati da spazi e 
stampa sullo standard output una riga di 23 numeri, di cui i primi 3 sono la ripetizione delle 
coordinate e gli altri sono le energie adiabatiche e diabatiche e le derivate parziali, come elencato 
in tab.3.2. 



e(l)- 


- X 


e(2)- 


- P3 


e(3)*- 


- z 


e(4)^ 


Eso 


e(5) ^ 


Es, 


e(6)^ 


dEso 


e(7)- 


\^ 



e(8)< 
e(9)< 
e(10) 
e(ll) 



dE^ 



So 



dPs 

dEs, 

dPs 

dEso 

dZ 

dEs, 



e(13) 
e(14) 



dZ 

e(12) ^ Vii (dial) 



e(15) ^ 

e(16) ^ 

e(17) - 

e(18) ^ 



V22 (dia2) 
V12 (coup) 



dx 
ay22 

dx 

ayi2 

dx 

dPs 



e(19) 
e(20) - 
e(21) 
e(22) 
e(23) 



ay22 

dPs 

aVl2 

dP -3 

dVn 
dZ 

ay22 

~dz' 

ayi2 

dz 



Tabella 3.2: Significato degli elementi del vettore e nella subroutine episa 

Commentando le linee comprese tra Begin gommano line client , End gommano line client il codice 
può essere compilato come oggetto esterno per importarne la subroutine episa mediante link. La 
subroutine richiede in input le coordinate X,P3,Z ed un vettore vuoto di dimensione 20, e, nel 
quale memorizzerà i valori elencati nelle posizioni in tab.3.2. 
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3.3 Funzione analitica per gli stati adiabatici 



Le variabili x,y,z corrispondono alle nostre X,P3,Z. 

! ############################# 
! Begin GOMMANO LINE CLIENT 
program pisa 

! E-Pisa - approximate energy v.3 
only for test purpose 
ci usage : 
$ . /episa 

$ echo "90 90 1" I ./episa 
real*8 x,y,z,t,a,e 
dimension e (20) 

13 format (f5.1,X,f5.1,X,f5.1,3X,2(fll.7,X)) 
15(f 10 .6,X)) 
e=0. dO 

write(6,*) "Episa command line" 
write(6,*) "Feed with Tors i on ( degrees ) , & 
Piramidalization (degrees ) , Stretching " 
write(6,*) "" Output in AU 
write(6,*) "x y z SO SI SO'x Sl'x SO'y & 
Sl'y SO'z Sl'z dial dia2 coup & 
dial 'x dia2 'x coup 'x dial 'y & 
dia2 'y coup 'y dial 'z dia2 ' z coup 'z" 
do while ( 1 . eq . 1) 
read (5 , *) x , y , z 

!conversion to internai coordinates 
z=z*0. 23477777 1387805 17 
cali episa (x , y , z , e) 
write(*,13) x,y,z,e 
enddo 
end 

End COMMAND LINE CLIENT 
############################# 
############################# 

MAIN SUBROUTINE 
############################# 

For external programs . 

Coupling of diabatic energies 

to obtain the adiabatic 

energies and the first 

derivatives in x,y,z. Feed the 

subroutine with real*8 <x,y,z> 

and <e> vector with dimension 20. 

Episa will fili <e> with these values : 

ADIABATIC VALUES: 

e(l)=SO, e(2)=Sl 

e(3) = S0'x, e(4)=Sl 'x 

e(5) = S0 'y, e(6)=Sl 'y 

e(7) = S0'z, e(8)=Sl 'z 

DIABATIC VALUES: 

e(9,10,ll)=dial ,dia2 , coup 

e (12 ,13 ,14)= (dial ,2, coup) 'x 

e (15 ,16 ,17)= (dial ,2, coup) 'y 

e(18,19,20)=(dial,2, coup) 'z 
subroutine episa (x , y , zi , e) 
real*8 x , y , z , zi , e , kcal2au 
dimension e (20) 

real*8 pi , vi , v2 , vl2 , autol , auto2 
real*8 vld , v2d , vl2d , autold , auto2d 
kcal2au=627.51d0 

!conversion to episa coordinate 
z=z 1/0. 23477777 1387805 17 
d2pi=dacos ( - 1 . dO ) /180d0 
vl=dial (x , y , z) 
v2=dia2 (x ,y ,z) 



vl2=coup (x , y , z) 
e(9)=vl 
e(10)=v2 
e(ll)=vl2 

autol=(vl+v2)/2 . dO 
auto2=((vl + v2)**2 
e(l)= autol - (auto 
e (2) = autol + ( aut 
! First Deriv 
vld=dialdx(x,y,z) 
v2d=dia2dx(x,y,z) 
vl2d=coupdx(x,y,z) 
e(12)=vld 
e(13)=v2d 
e(14)=vl2d 
autold=(vld+v2d)/2 . 
auto2d=0.25d0*(auto 
auto2d=auto2d*(2 . dO 
4. d0*(vl*v2d+vld*v2 
e (3) = autold - auto2d 
e (4) =autold+auto2d 
! First Deriv 
vld=dialdy(x,y,z) 
v2d=dia2dy (x,y ,z) 
vl2d= coupdy (x , y , z) 
e(15)=vld 
e(16)=v2d 
e(17)=vl2d 
autold = (vld + v2d)/2 . 
auto2d=0.25d0*(auto 
auto2d=auto2d*(2 . *( 
4. *( vl*v2d+vld*v2) 
e (5) = autold - auto2d 
e (6)=autold+auto2d 
! First Deriv 
vld=dialdz(x,y,z) 
v2d=dia2dz (x,y ,z) 
vl2d= coupdz (x , y , z) 
e(18)=vld 
e(19)=v2d 
e(20)=vl2d 
autold=(vld + v2d)/2 . 
auto2d=0 . 2 5d0*(auto 
auto2d = auto2d*(2 . dO 
4. *( vl*v2d+vld*v2) 
e (7) = autold - auto2d 
e (8)=autold+auto2d 
do i=l,20 
e(i)=e(i)/kcal2au 
enddo 

e(7)=e(7)/0. 2347777 
e(8)=e(8)/0. 2347777 
e(18)=e (18)/0. 23477 
e(19)=e (19)/0. 23477 
e(20)=e(20)/0. 23477 
end 



) + 4. d0*(vl2**2) -4. d0*vl*v2 
2**0.5d0)/2. dO; 
o2**0.5d0)/2. dO; 
ative in X 



dO 

2)**(-0.5d0) 
*(vl+v2)*(vld+v2d)- 
) + 8. d0*vl2*vl2d) 



ative in Y 



2)**(-0.5d0) 
vl+v2)*(vld+v2d) ■ 
+ 8. d0*vl2*vl2d) 



ative in Z 



2)**(-0.5) 
*(vl+v2)*(vld+v2d) - 
+ 8. d0*vl2*vl2d) 



7138780517d0 

7138780517d0 

777138780517d0 

777138780517d0 

777138780517d0 



Seguono in totale 12 definizioni di funzioni nella forma: 

! Diabatic surfaces 

real*8 function dial(x,y,z) 

real*8 x, y, z 

dial= . . . 

return 

end 



Tabella 3.3: il programma Ep 
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4 Dinamica semiclassica 



Per investigare il meccanismo di fotoisomerizzazione del CPP, analizzare l'effetto del moto delle 
singole coordinate nucleari sulla sua resa quantica ed infine per determinare il modo in cui la 
reazione ha luogo dopo la fotoeccitazione, è possibile condurre simulazioni di tipo dinamico. 

Studi di questo tipo si possono svolgere secondo diversi approcci e con diverse accuratezze. 
L'approccio meno dispendioso computazionalmente è quello della dinamica molecolare classica, 
che per questa caratteristica è il metodo di adozione per studi su molecole di grandi dimensioni, 
come proteine. 

Il moto classico della molecola è confinato su una singola superficie elettronica adiabatica (quella 
fondamentale), impendendo la descrizione di transizioni non radiative innescate da accoppiamenti 
non adiabatici. 

Questo approccio non può essere applicato nel nostro caso, poiché la reazione che intendiamo 
studiare inizia sullo stato eccitato Si per poi procedere sullo stato fondamentale So in seguito ad 
una transizione non radiativa. 

All'altro estremo del piano costo-accuratezza vi è la dinamica quantistica, metodo dispendioso 
computazionalmente ma molto preciso nei risultati. 

Il metodo che abbiamo scelto per il nostro studio è quello semiclassico, utilizzato ad esempio 
da Ferretti et al. [10], ovvero un ibrido tra dinamica classica e quantistica. La molecola si muove 
prevalentemente in modo classico Newtoniano sulle superfici di energia potenziale; l'evoluzione 
dell'equazione di Shrodingrer dipendente dal tempo è considerata con la tecnica di Lanczos [12, 
15,30,31] ed il passaggio da uno stato elettronico all'altro è governato dall'algoritmo di surface 
hopping di Tully, denominato fewest switches [51,52], brevemente esposto di seguito. 

Il metodo fewest switches Consideriamo due stati elettronici e 1. Lo stato dipendente dal 
tempo è così definito: 

m)) = Mt)e-'^^'\MQ) (4.1) 

Dove Q è il vettore delle coordinate cartesiane dei nuclei e la fase 70 è data da: 

lo{t) = \j EoiQit'))dt' 
Sostituendo 4.1 nell'equazione di Shcrodinger: 

in^^ = H,,mt)) (4.2) 
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4.1 Definizione del pacchetto d^onda iniziale 
si ottiene un sistema di equazioni differenziali per l'ampiezza di probabilità nello stato 0, Aq: 

Ào{t) = -A,{t)e^^^-^Y.Q»^^<y>l (4.3) 

a 

7o,i =7o -7i 

Gq[ è una componente del vettore di accoppiamento non-adiabatico calcolato dalla derivata 
delle energie adiabatiche H rispetto dalle coordinate nucleari Qa, proiettata sulla base diabatica 

ÌVaìAVb)- 

m = (1%), |%))Co IV'i) = (1%), |%))Ci (4.4) 



Gt\ = {i^o,i 



d 



dQa 



i;,) = (Eo-E,)-'cl^C^ (4.5) 



Le probabilità adiabatiche semiclassiche si ottengono dal modulo quadro dell'ampiezza di pro- 
babilità Po{t) = |Ao(t)p. L'algoritmo fewest switches di Tully stabilisce che, in un problema a 
due stati, la probabilità di salto da ìpo ^ i^i è data da: 



Wo^i = ^At, se Po > 0, 

V|/o,i = 0, se Po <= (4.6) 

Un numero casuale stabilisce se il salto debba avvenire o meno, a seconda che sia maggiore 
o minore della probabilità VFo,i- Nel caso in cui un salto abbia luogo l'energia cinetica viene 
modificata per assicurare la conservazione dell'energia, modificando i momenti lungo la direzione 
del derivative coupling in eq.4.5. 

Il metodo è stato implementato in linguaggio Fortran da F. Santoro e adattato nel corso del 
lavoro di tesi per includere il nostro potenziale (sez.3.3.1). Ci riferiremo a questa implementazione 
col nome TSH . 



4.1 Definizione del pacchetto d'onda iniziale 

L'assorbimento di un fotone da parte del CFP porta la molecola dallo stato elettronico fonda- 
mentale a quello eccitato mantenendo cioè la geometria nucleare e lo stato vibrazionale propri 
dello stato di partenza, in accordo con il principio di Franck Condon. La posizione ed i momen- 
ti nel punto Frank-Condon rappresentano dunque il punto di partenza per qualsiasi dinamica: 
da lì la geometria molecolare evolverà sulla superficie dello stato eccitato lungo una traiettoria 
prevedibile, fino al decadimento. 

L'impulso luminoso, per quanto breve, eccita un gruppo di molecole aventi ciascuna una po- 
sizione ed un momento entro una certa distribuzione statistica attorno al punto di equilibrio, 
formando così un pacchetto d'onda. Per valutare il comportamento di un grande numero di 
molecole abbiamo definito un pacchetto d'onda iniziale composto da 10^ molecole, con posizioni 
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4 Dinamica semiclassica 

e momenti iniziali calcolati secondo la distribuzione di Wigner del modulo quadro della funzione 
d'onda dello stato fondamentale, nello spazio delle coordinate XPZ e dei relativi momenti. 

La funzione d'onda impiegata è una approssimazione armonica del potenziale nel punto di 
minimo Trans: le frequenze di oscillazione armonica attorno al punto di equilibrio derivano dal 
fit dell'energia potenziale con una parabola (4.1). 



Coordinata: 



Frequenza Armonica: 



X 



69.24 cm- 



341.42 cm- 



1332.75 cm- 



Tabella 4.1: Frequenze di vibrazione lungo le tre coordinate 

La distribuzione di posizioni e momenti ottenuta segue l'andamento di una Gaussiana (fig.4.1), 
sia nello spazio delle coordinate che in quello dei momenti. 






■• •''.9vj:^*;.j:«:'>'7»'V.-. ' 



Figura 4.1: Distribuzioni di Wigner di "^frans P^^ ^^ posizioni iniziali 

L'evoluzione della molecola viene calcolata a partire da ciascun punto del pacchetto separata- 
mente. Poiché esso è molto localizzato attorno alla geometria dello stato fondamentale, il pac- 
chetto d'onda si muoverà sulle PES inizialemente come un oggetto fisico coerente. Nel corso del 
tempo la piccola distribuzione del pacchetto iniziale viene amplificata, ed il pacchetto si disperde 
sempre più lungo le tre coordinate, fino a non poter essere trattato come un oggetto coerente. 

4.2 Verifica del fewest switches 

Il calcolo è stato svolto su 10^ traiettorie, senza attrito, con una evoluzione di 2ps divisa in step 
di integrazione delle equazioni differenziali del moto di 0.5fs. 

L'algoritmo di fewest switches calcola la probabilità che una traiettoria salti da uno stato all'al- 
tro in modo che, per un numero di traiettorie statisticamente significative, le frazioni di traiettoria 
sui due stati corrispondano alle probabilità adiabatiche derivanti dalla soluzione dell'equazione di 
Scrodinger dipendente dal tempo [33]. In fig.4.2 sono confrontate le probabilità con le frazioni di 
traiettorie sui due stati ad ogni tempo. L'accordo tra i due dati è ottimo, e ciò conferma che il 
fewest switches viene correttamente applicato. 

Il calcolo ha tuttavia presentato un numero di salti tra uno stato e l'altro molto elevato: la 
maggior parte delle traiettorie cambia stato tra le 6 e le 10 volte in 2ps. Se analizziamo il valore 
di torsione al quale è avvenuto l'ultimo salto (US) prima del termine della traiettoria, troviamo 
una distribuzione molto ampia lungo i valori di torsione 4.3(a). Neppure il 25% delle traiettorie 
ha l'US entro 10° dall'intersezione, mentre 15% ha l'US entro 10° dalle geometrie Trans e Cis 
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4.2 Verifìca del fewest switches 



Contronto tra probabilità e frazioni di 
traiettorie 




le+03 1.25e+03 

Tempo (fs) 



r igura 4.2: Confronto tra probabilità adiabatica e frazioni di traiettorie sui due stati. 

dove i due stati sono molto lontani e la probabilità di transizione dovrebbe essere, fisicamente, 
molto ridotta. Il resto degli US è distribuito lungo tutti gli altri valori di torsione. 




(a) Numero di salti US, sommati lungo la torsione tra X ed (b) Numero di salti US, sommati lungo la torsione tra X ed 
X+l° X+l°, dopo lo smooth sull'accoppiamento 

Figura 4.0! Risoluzione dell'anomalia nel numero di salti in funzione della torsione. 



Questa anomalia si spiega considerando che la molecola trascorre molto tempo sui punti Trans 
e Cis rispetto a quello trascorso in prossimità dell'intersezione conica (vedi sez.4.4.1). Seppure 
gli accoppiamenti non-adiabatici GqJ^ siano molto piccoli e così VFo,i, le traiettorie hanno molto 
tempo per saltare tra i due stati. L'algoritmo fewest switches viene interrogato con numeri casuali 
molte più volte presso Trans e Cis, e questo bilancia in qualche modo il piccolo VFq,! portando un 
risultato poco fisico. 

In TSH abbiamo perciò applicato uno smooth Gaussiano agli accoppiamenti non-adiabatici: 






)/{50kcal/mol))^ 



(4.7) 



In tal modo le probabilità Po{t) sono rapidamente soppresse in regioni dove gli stati adiabatici 
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4 Dinamica semiclassica 

sono troppo diversi, e la distribuzione di salti risultante 4.3(b) è quella fisicamente ragionevole, 
con la maggior parte delle transizioni localizzate attorno all'intersezione conica. 

Questo ha un effetto negativo sulla coerenza tra le frazioni di traiettorie sugli stati, frutto del 
fewest switches, e la probabilità adiabatica derivante dall'equazione di Schrodinger. Dalla fig.4.4 
si rileva infatti una discrepanza di circa .2 tra i due dati, che abbiamo trascurato in favore della 
maggior fisicità dei salti. 

Confronto probabilità, frazioni dopo lo 
snnoothing 
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Tempo (fs) 



Figura 4.4: Confronto tra probabilità adiabatica e frazioni di traiettorie sui due stati dopo l'applicazione dello smoothing 



4.3 Dinamica con smoothing e stop 

Introdotto lo smoothing in eq.4.7 abbiamo condotto un nuovo calcolo su 10"^ traiettorie, con tempo 
finale di 2ps, e step di integrazione di .5fs. 

Su ogni traiettoria è stata imposta una condizione di uscita sul valore di torsione allo stato 
fondamentale. Ad ogni step di integrazione il programma controlla se la traiettoria si trova sullo 
stato fondamentale, quindi, se | cosX| < cosS^"'^ la traiettoria viene fermata ed il calcolo passa a 
quella successiva. Questo comportamento simula in modo molto semplice e computazionalmente 
efficiente la presenza dell'attrito, che nel corso del tempo fermerebbe i pacchetti nei minimi dello 
stato fondamentale. 

Ad ogni step di integrazione TSH abbiamo riportato le popolazioni sui due stati adiabatici, 
classificate come Trans per 90°<X<270° e Cis per X>270°, X<90°. Al fine di ottenere un insieme 
di dati più omogenei e comprensibili riguardo la coordinata X, cui il potenziale periodico permette 
di assumere qualsiasi valore positivo o negativo, abbiamo proiettato nell'intervallo 0-360° il valore 
calcolato dalla soluzione del moto. 

La figura 4.5 mostra come fino a 400fs non ci sia sostanziale decadimento dallo stato eccitato 



""^l cosX| < cos5° è verificata ogni volta che X assume, in modulo 360, valori -5<X<5, 175<X<185, ovvero intorno 
ai minimi lungo la torsione 
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4.3 Dinamica con smoothing e stop 
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Figura 4.5: Popolazioni adiabatiche Cis e Trans sui due stati 

a quello fondamentale, né alcuna oscillazione tra conformazioni Trans e Cis su Si. Dopo questo 
lungo lasso di tempo dove il sistema è apparentemente immobile, nel giro di 200fs più del 70% 
delle traiettorie lascia la geometria Trans. Di queste circa la metà decade allo stato fondamentale 
So, e l'altra metà prosegue invece sullo stato eccitato per assumere conformazione Cis. 

Al termine della dinamica circa il 90% delle molecole è decaduto e di queste il 60% ha assunto 
geometria Cis, quasi il 30% si è attestato invece a Trans. Ridistribuendo il restante 10% di 
traiettorie ancora sullo stato eccitato secondo le stesse proporzioni si ottiene il 66% di molecole 
Cis ed il 33% di Trans. In altre parole la resa quantica è del 66%, in linea con quanto misurato 
per la Rodopsina [35]. 

In fig.4.6 sono invece riportate le popolazioni sugli stati diabatici in funzione del tempo. Le 
rapide oscillazioni sono dovute al moto lungo lo stretching, come si vede in fig.4.7(d),4.7(c): la 
coordinata Z oscilla sulla superficie di energia potenziale dello stato adiabatico, attraversando i 
due stati diabatici ad ogni periodo. 

Nelle figure 4.7 è riportata la media delle coordinate P3 e Z ad ogni step di integrazione su 
tutte le traiettorie. La coordinata X non è mostrata poiché la simmetria delle superfici ne rende 
non significativo il dato. 

La piramidalizzazione media si comporta in modo simile su entrambi gli stati, seguendo il 
cammino di energia minima già individuato in fig.3.2(c). Su Si la coordinata P3 diminuisce di 
circa 1° nel corso dei primi 400fs, quindi torna a crescere oltrepassando il valore di equilibrio. 
Dopo Ips il pacchetto d'onda si disperde, restituendo un valor medio instabile numericamente. 
Su So il numero di traiettorie non è statisticamente significativo fino a 400fs, per quanto visto in 
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4 Dinamica semiclassica 



Popolazioni diabatiche 
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Figura 4.6: Popolazioni sui due stati diabatici 

fig.4.5; successivamente la coordinata rimane stabile al valore di equilibrio. 

Lo stretching su Si presenta una sorprendente coerenza, tanto che a 2ps il valor medio oscilla 
ancora molto nitidamente. L'ampiezza di oscillazione comincia a decrescere lOOfs prima del - ed 
in concomitanza con - il trasferimento di popolazione allo stato fondamentale: questo suggerisce 
che inizialmente l'energia di cinetica sia in gran parte immagazzinata nello stretching, il quale 
comincia a cederla agli altri gradi di libertà dopo 200fs. 

4.4 Pacchetti d'onda e traiettorie sullo stato eccitato 

Durante il calcolo è stata seguita l'evoluzione dei pacchetti d'onda memorizzando la posizio- 
ne di tutti i pacchetti ogni 50fs; di 100 traiettorie si è tracciata la posizione ogni 4fs fino al 
raggiungimento dei 2ps o allo stop sullo stato fondamentale. 

Nelle fig.4.8 sono disegnati i pacchetti come punti in funzione di due coordinate per volta. In 
fig.4.8(c),4.8(b) si nota il pacchetto altamente localizzato su in tre valori di stretching diverso a 
50,150,200fs (rosso, arancione, giallo): come suggerito dalle coordinate medie, il pacchetto si muove 
coerentemente lungo Z. Si confronti in particolare la quasi totale dispersione raggiunta lungo la 
torsione in tempi lunghi (punti viola e neri) in fig.4.8(c) con la distribuzione lungo Z ancora 
compatta e ben riconoscibile. 

Lungo la torsione vediamo che le il pacchetto rimane almeno fino a 300fs localizzato attorno 
al punto di Franck Condon di 180° senza disperdersi significativamente (rosso, arancione, giallo, 
verde e azzurro): rimane sostanzialmente immobile lungo X almeno fino a 300fs (verde). 
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4.4 Pacchetti d^onda e traiettorie sullo stato eccitato 




Tempo (fs) 
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(a) Piramidalizzazione P3 inedia sullo stato eccitato Si 



Piramidalizzazione media su SO 
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(b) Piramidalizzazione P3 media sullo stato fondamentale 

So 



Stretching medio su 51 
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(c) Stretching medio Z sullo stato eccitato Si 



stretching medio su SO 
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(d) Stretching medio Z sullo stato fondamentale So 



Figura 4.7: Coordinate mediate su tutte le 10^ traiettorie, in funzione del tempo. Stato eccitato Si a sinistra, fondamentale 
So a destra. 
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4 Dinamica semiclassica 

Lungo P3 il pacchetto si dispone inizialmente intorno e sotto i 90° e sotto (vedi verde e azzurro), 
quindi dopo 450fs raggiunge anche piramidalizzazioni superiori a 90°. 



(a) Torsione in ascissa e piramidalizzazione in ordinata. 



r^ii4^]^^ife: i '■' 
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(b) P3 in ascissa e Z in ordinata. 
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(e) Torsione in ascissa e stretching in ordinata. 

Figura 4.8: Pacchetto d'onda a: rosso=50fs; arancione= 150; giallo= 200; verde= 300; azzurro= 450; viola= 750; blu= 
1000; grigio= 1250; nero= 1650fs 

Le 100 traiettorie campionate ogni 4fs contribuiscono a chiarire il quadro dinamico, assieme 
agli integrali di tutti i pacchetti entro un certo intervallo, fig.4.9. 

La coerenza del moto lungo Z (fig.4.9(e)) è evidenziata dall'andamento delle 100 traiettorie 
sovrapposte, che su questa coordinata si comportano quasi fossero una sola per i primi 600fs 
circa. Riguardo la piramidalizzazione non otteniamo informazioni aggiuntive dal grafico 4. 9 (e), 
se non una conferma della decoerenza del moto lungo P3. 

Da fig.4.9(a),4.9(b) emerge che la torsione rimane effettivamente immobile per i primi 400fs: 
durante questo lungo periodo di induzione continua a muoversi regolarmente sulle altre due 
coordinate, in particolare lungo lo stretching tra -.5 e 3. 
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4.4 Pacchetti d^onda e traiettorie sullo stato eccitato 




(a) Torsione X in funzione del tempo 



(b) Numero di pacchetti per intervallo di torsione (1°) 





(e) Piramidalizzazione P3 in funzione del tempo (d) Numero di pacchetti per intervallo di piramidalizza- 

zione (1°) 





(e) Stretching Z in funzione del tempo 



(f) Numero di pacchetti per intervallo di stretching (.1) 



Figura 4.9: Le coordinate X,P3,Z del pacchetto d'onda; a sinistra, monodimensionale con dati campionati ogni 4fs; a 
destra grafico 3D del numero di pacchetti entro un certo intervallo della coordinata (campioni ogni 50fs; su 

tutte le traiettorie). Nota: le traiettorie che vanno a che compaiono da repentinamente hanno saltato sullo/dallo stato 



fondamentale 
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4 Dinamica semiclassica 

4.4.1 Trapping delle traiettorie 

L'esteso tempo di induzione necessario per muovere la torsione dal punto di Franck Condon ed 
avviare l'isomerizzazione è stato analizzato da Olivucci, Santoro e Lami in [38,43]. L'effetto 
è dovuto alla presenza di un plateau (fig.4.10(e)) o addirittura di un minimo lungo la torsione 
(fig.4.10(c)) nella geometria Franck Condon sullo stato eccitato. 

Dopo l'eccitazione la molecola non ha alcuna accelerazione torsionale (fig.4.10(d)), mentre la 
PES ha una forte pendenza verso valori di stretching più alti (1.5). Inizialmente dunque la 
molecola accelera lo stretching, fino ad arrivare in regioni dove il plateau o il minimo scompare, 
sostituito da un massimo (fig.4.10(f)). A questo punto il potenziale comincia a muovere lentamente 
la torsione (la frequenza di vibrazione è molto bassa) ma tutta l'energia cinetica è quasi tutta 
sullo stretching, che continua il suo moto verso valori che eccedono il punto di equilibrio di Z. 

Intorno a Z=2.7 ricompare (fig.4.10(f)) un minimo locale a torsione 180°, che diviene pia profon- 
do all'aumentare di Z (fig.4.10(h)). Il moto di stretching si ripete quindi in senso opposto: ad ogni 
oscillazione dello stretching, il minimo locale lungo la torsione la focalizza a 180°, intrappolando 
così la molecola sullo stato eccitato. 

Esistono diverse strategie per affrontare questo fenomeno di trapping: 

• Popolare gli stati termici ad energia più elevata, in modo da rendere cineticamente insigni- 
ficante la piccola barriera energetica che intrappola la molecola. 

• Estendere la distribuzione iniziale dei momenti, con conseguenze analoghe (più energia nel 
sistema) 

• Svolgere un calcolo di dinamica quantistica. 

In [43] gli autori osservano un comportamento caotico per calcoli di dinamica classica su PES con 
larghi plateau o piccoli minimi: il risultato della dinamica dipende enormemente dalle condizioni 
iniziali che si scelgono per il pacchetto. La propagazione quantistica del pacchetto d'onda risolve 
invece gli artefatti introdotti dal plateau, permettendo di stimare il reale effetto dei minimi locali 
sul decadimento. 
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4.4 Pacchetti d^onda e traiettorie sullo stato eccitato 
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r igura 4.10: Effetto trapping torsionale per valori di stretching molto minori del punto di equilibrio 
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4 Dinamica semiclassica 

4.5 Conferma al modello ITP 

Al termine di ogni traiettoria sullo stato fondamentale è memorizzata la geometria finale assieme 
alle coordinate a cui è avvenuto l'ultimo salto dallo stato eccitato, detto salto reattivo. 

Abbiamo così modo di verificare mediante la dinamica semiclassica il modello di interazione 
torsione-piramidalizzazione introdotto qualitativamente in sez.2.5. La figura 4. 11 (e) mostra che 
la maggior parte dei salti reattivi verso Trans avviene a piramidalizzazione P3 maggiore di 90° di 
5-10°, in particolare attorno all'intersezione conica (X=90°,270°). 

In figura 4. 11 (e) la maggior concentrazione di salti reativi a Cis si trova invece, presso la CI, a 
piramidalizzazioni inferiori a 90° di 5-10°. 

Nelle figure 4.11(f),4.11(d) il numero di salti reattivi che portano a Cis/Trans è integrato su 
intervalli di piramidalizzazione di 1° e stretching di 0.1. I salti reattivi con decadimento Cis 
sono nettamente concentrati su piramidalizzazioni inferiori a 90° (come previsto dal modello 
ITP); quelli con decadimento Trans sono più regolarmente distribuiti attorno a 90°, con una coda 
dissimetrica verso valori superiori a 90°. 

La correlazione con Z non è così netta tra decadimento Cis e Trans, anche se è osservabile una 
maggior concentrazione di salti reattivi Cis per Z<1 e viceversa per Z>1. 
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2 Z 



(a) I punti indicano la posizione dell'ultimo salto Si -^ Sq che ogni tra- 
iettoria ha eseguito prima del termine della simulazione. In rosso sono 
colorati quelli che si trovano in regione Cis, in blu quelli in regione Trans. 






(e) Salti con decadimento Trans. (d) Integrazione sullo stretching dei salti reattivi. In Blu i salti reattivi 
Torsione in ascissa, piramidalizza- che riportano a Trans; in rosso quelli che conducono a Cis. 
zione in ordinata. 

Ntra 
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(e) Salti con decadimento Cis. Tor- (f) Integrazione sulla piramidalizzazione dei salti reattivi. In Blu i 

sione in ascissa, piramidalizzazione salti reattivi che riportano a Trans; in rosso quelli che conducono a 

in ordinata. Cis. 

Figura 4.11: Correlazione ultimo salto-geometria di termine 
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5 Conclusioni 



Nel corso di questo lavoro di tesi abbiamo: 

• Calcolato le superfici di energia potenziale sugli stati elettronici fondamentale So ed eccita- 
to Si per il 4-(cyclo-pent-2-enyldiene)-3,4-dihydro-2H-pyrrolinium (CPP) in un insieme di 
coordinate rilevanti per la fotoisomerizzazione (C^ea), composte da una coordinata torsio- 
nale, una di piramidalizzazione sull'anello pirrolinico ed una di stretching alternato lungo il 
sistema tt. 

• Individuato una interazione tra le coordinate di torsione e piramidalizzazione (ITP), de- 
terminante linee di quasi degenerazione tra i due stati elettronici e regioni di rilassamento 
Cis/Trans sullo stato fondamentale. Questa interazione ha conseguenze dirette sulla resa 
quantica di reazione, confermate da simulazioni di dinamica molecolare semiclassica. La 
causa dell'ITP è stata riconosciuta nella morfologia degli orbitali molecolari dello spazio 
attivo definito per il calcolo CASSCF. 

• Diabatizzato le superfici adiabatiche utilizzando il momento di dipolo elettrico sui due stati 
elettronici come indice dello stato diabatico molecolare. Le superfici diabatiche risultanti 
e l'accoppiamento tra di esse sono smooth rispetto alle coordinate nucleari, permettendo 
calcoli di dinamica più accurati. La diabatizzazione realizzata ha senso fisico maggiore 
rispetto a quelle ottenute con altri metodi proposti in letteratura per sistemi simili. 

• Condotto una prima simulazione di dinamica molecolare semiclassica, la quale ha permesso 
di identificare una caratteristica topologica della PES sullo stato eccitato che promuove il 
trapping di traiettorie sullo stato eccitato per tempi molto lunghi. I risultati inoltre hanno 
confermato il modello di interazione torsione-piramidalizzazione emerso dall'analisi delle 
PES e degli orbitali molecolari. 

Il raggiungimento di questi obiettivi ha comportato l'acquisizione di conoscenze di meccanica 
molecolare, come la definizione di coordinate interne e la compilazione di una OJtairice Z , e 
di calcolo elettronico CASSCF con i programmi Qaussian e -^lcas . L'estrazione dei dati dai 
file di output, la manipolazione di OJlairice Z e la distribuzione del calcolo su cluster virtuali 
attraverso rete TCP/IP, anziché essere considerati compiti gravosi e ripetitivi, sono stati affrontati 
ricercando la massima automazione, e colti come opportunità per approfondire la padronanza del 
linguaggio di programmazione Python. L'elaborazione della mole di dati ottenuta mi ha consentito 
di acquisire dimestichezza con il programma di algebra computazionale Mathematica. 
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5.1 Sviluppi futuri 

5.1 Sviluppi futuri 

Il lavoro di ricerca futuro consisterà nello svolgimento della dinamica quantistica, grazie all'aiuto 
di Alessandro Lami, per confrontarne i risultati con la dinamica semiclassica. 

L'algoritmo di switch per la dinamica semiclassica dovrà quindi essere semplificato senza alterare 
le proprietà dinamiche del sistema, al fine di implementarne una versione meno esigente di risorse 
computazionali nel modello di proteina sviluppato dai nostri collaboratori all'EMPA [6]. 

Il modello di interazione torsione-piramidalizzazione è promettente per il design in silico di nuovi 
switch molecolari. L'estensione dell'ITP a tutti i gradi di libertà ne verificherà gli effetti sulla resa 
quantica di reazione, aprendo nuove prospettive di ricerca nella comprensione dell'efficienza dei 
processi fotochimici e nella sintesi di nuove specie fotoattive. 
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A ReQuieM 



L'acronimo fuzzy^ ReQuieM sta per Remote Queue Manager, ovvero Gestore di Code 
Remote. Il programma è nato per rispondere all'esigenza di condurre una scansione Qaussian in 
uno spazio di tre o più dimensioni e con il vincolo, esposto in sez.2.3, di lanciare il calcolo in un 
punto di questo spazio solo se si posseggono già gli orbitali per il CASSCF di un punto attiguo. 
ReQuieM è stato poi esteso per consentire la distribuzione del calcolo su più computer e 
cluster di computer. L'Istituto per i Processi Chimico-Fisici di Pisa ha infatti reso disponibili per 
questo lavoro di tesi 4 macchine biprocessore ( cui ci riferiremo con i nomi mir, ciro, toto, rosa) 
e 6 processori divisi in due account sul cluster (hydra). Era prevista la possibilità di accedere ad 
un elevato numero di processori su un cluster presso il Center for Photochemical Sciences della 
Bowling Green State University; infine Olivucci aveva prospettato la disponibilità di ulteriori CPU 
su cluster presso il laboratorio di calcolo di Siena ed il CICAIA di Bologna. 

La potenza di calcolo ottenibile da tutti questi processori avrebbe consentito di mappare l'e- 
nergia potenziale del CPP con una grande risoluzione, tuttavia si trovavano su computer e cluster 
di computer differenti, ciascuno con propri codici di accesso, sistema di gestione, disponibilità 
di calcolo. Lanciare e seguire manualmente i calcoli avrebbe richiesto un operatore pronto ad 
intervenire 24h/24h, pena l'allungamento smisurato dei tempi di scansione. In ReQuieM sono 
allora state inclusi metodi per inviare lavori presso account remoti, controllarne l'avanzamento ed 
infine recuperarne i risultati, senza alcun intervento manuale. 

Il programma opera su un nodo master (o locale) che necessita di poca capacità di calcolo e 
molta di immagazzinamento dati. Lo schema di esecuzione è il seguente: 

I Stabilisce una connessione criptata con tutte le macchine configurate. 

II Organizza le scansioni da effettuare in spazi di scansione^ in strutture dati opportune. 

Ili Itera le seguenti operazioni fino al termine della scansione: 

1 Legge un file di comunicazione per aggiungere o rimuovere spazi di scansione. 

2 Calcola e scrive alcune statistiche circa l'avanzamento della scansione. 

3 Ricerca nello spazio di scansione un punto non completato che soddisfi la condizione 
di attiguità con almeno un punto completato. 

4 Ricerca un processore libero tra tutti i processori configurati. 



^11 singolare acronimo impreciso nasce durante l'ascolto dell'omonima opera di Mozart per analogia tra la stupe- 
facente polifonia iniziale dell'Introitus e le molteplici risorse computazionali allocate dal programma, a partire 
dal job iniziale per propagarsi esponenzialmente a quelli attigui, come la massa sonora si muove dai bassi agli 
alti nel primo tema della partitura. 
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5 Avvia il calcolo del punto sul processore libero, inviando alla macchina su cui risiede 
tutti i file necessari. 

6 Controlla lo stato di tutti i calcoli in esecuzione. 

7 Recupera i risultati dei calcoli conclusi, trasferendoli dalle rispettive macchine al nodo 
locale. 

8 Salva tutte le informazioni riguardanti il proprio stato (snapshot). In tal modo è 
possibile interrompere e riprendere la scansione. 

9 Attende, nel caso in cui tutte queste operazioni non abbiano richiesto almeno un certo 
tempo. 

IV Salva un'ultima volta il proprio stato; chiude le connessioni remote; scrive le statistiche di 
scansione; saluta l'operatore ed esce. 

ReQuieM è implementato in Python [39], un linguaggio di programmazione dinamico orien- 
tato ad oggetti molto versatile, con regole sintattiche che migliorano leggibilità, manutenibilità e 
riusabilità del codice. È dotato di ottima documentazione e di una ricca libreria standard; mol- 
te estensioni sono reperibili facilmente in rete. In particolare la connettività SSH [50] è fornita 
dal pacchetto esterno Pexpect [37], che rende le comunicazioni criptate semplici come leggere e 
scrivere su un file. L'utilissima e indispensabile funzionalità di snapshot deriva invece dal modulo 
cPickle della libreria standard di Python. 

Il codice di ReQuieM è suddiviso in quattro moduli, rappresentati da altrettanti file Python 
.py e chiamati con il nome senza estensione: glassian , requiem , ClusterConf e RequiemDirector 
La modularizzazione semplifica notevolmente: 

• la manutenzione del codice: ogni modulo contiene funzionalità specifiche e rimane limitato 
ad alcune centinaia di righe 

• il debug, essendo il nome del modulo specificato nelle eccezioni e nei traceback a seguito di 
crash 

Nelle prossime sezioni sono presentati i moduli con le classi di oggetti e/o alcune funzioni di 
particolare importanza per comprendere il funzionamento interno di ReQuieM . 

Formato della griglia 

Le griglie di punti cui ci siamo estesamente riferiti nei capitoli precedenti sono rappresentate 
nella memoria fissa del computer da cartelle di file, contenenti un file di input di Qaussian per 
ciascun punto. Una cartella deve avere un punto iniziale completo di output e checkpoint (file 
dove Qaussian memorizza gli orbitali), e può contenere solo punti collegabili al punto iniziale 
percorrendo solo punti attigui tra di loro: per questi punti infatti non è possibile soddisfare la 
condizione di attiguità e sono pertanto ignorati da ReQuieM . 

ReQuieM stabilisce a quale punto della griglia i punti appartegnono in base al nome che hanno 
sul filesystem i file di input che li rappresentano. I formato è xN-yM-zK.com, dove N,M,K sono 



85 



A ReQuieM 

numeri naturali che indicano lo step sul quale si trova il punto lungo la coordinata relativa (x,y,z). 
ReQuieM utilizza questi indici per stabilire se due punti a,b sono attigui: 

l^a - ^ò| + \Ma - M5I + \Ka - i^^l = 1 condizione di attiguità (A.l) 

Un file scan.info contiene la traduzione degli indici N,M,K in unità fisiche per permettere di 
raccogliere e presentare i risultati al termine della scansione. 

Griglie input Qaussian in questo formato sono facilmente generabili con GeoM<n . 

A.l Metadati, analisi dell'output: glassian 

glassian ^ contiene due classi ed alcune semplici funzioni per la manipolazione di stringhe. 

A. 1.1 class GaussianOutput: 

è la rappresentazione di un file di log di Qaussian e contiene metodi per il parsing e l'estrazione 
dei risultati di interessse, che sono tradotti in strutture dati Python. Contiene le funzioni: 

decode_name() decodifica il nome sul filesystem dell'input, ricavando gli indici N,M,K per le 
tre coordinate. 

duration() legge la durata del calcolo dal log 

energies() estrae le energie su tutti gli stati elettronici 

mulliken() legge le cariche di Mulliken sugli atomi e somma quelle appartenenti ad un anello o 
all'altro 

readlog(log) apre il file log ed estrae durata, energie, cariche di mulliken con i metodi precedenti. 

zmatrix() estrae la OJlairice Z del punto e la restituisce come stringa 

Gli oggetti sono inizializzati indicando un log da analizzare: obj=GaussianOutput ( Vpath/to/f ile.log') 

A. 1.2 class ScanLPS: 

La classe è stata solo impostata per ospitare la logica che determina il tipo di scansione tra quelli discussi 
in sez.2.3, ma per il momento è sostanzialmente vuota. Il metodo scaninfo() legge il file scan.info e assegna 
gli attributi xinfo, yinfo, zinfo, contenenti la traduzione degli indici dei punti in coordinate fisiche. 

A. 2 Opzioni di configurazione: ClusterConf 

ClusterConf contiene opzioni di configurazione sotto forma di assegnazioni di varibili e inizializ- 
zazioni di oggetti Machine (A. 3. 3). 



^11 nome di questo modulo è un gioco di parole glass (occhiale, lente) + Qaussian . Come una lente aiuta a 
decifrare scritture molto fìtte, glassian contiene la logica per eseguire il parsing dell'output assai verboso di 
Qaussian 
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A. 2 Opzioni di configurazione: ClusterConf 



1 maiii_logf ile= ' /home/ daniele /data/requiem, log ' 

2 main_loglevel=0 

3 maiii_logname= ' main ' 

4 main_iiiput_f ile= ' /home/ daniele /data/ input _ file ' 

5 main_dump= ' /home/ daniele /data/requiem . obj ' 

6 sleepness=5 
7 

8 f abh=Machine (name= ' fabh ' , 

9 maiiager=' queue ' , 

10 host= ' hydra ' , 

11 user='XXXXXXXX' , 

12 password='XXXXXXXX' , 

13 workload = 5, 

14 workdir= ' /home/ fabriz io / daniele /requiem/ ' , 

15 poll_interval=sleepness /5) 

16 lamh=Machine (name= ' lamh ' , 

17 manager= ' queue ' , 

18 host= ' hydra ' , 

19 user='YYYYYYYYY' , 

20 password='YYYYYYY' , 

21 workload = 5, 

22 workdir= ' /home/lami/ daniele /requiem/ ' , poll_interval=sleepness /5) 

23 ciro=Machine (name= ' ciro ' , 

24 manager= ' direct ' , 

25 host= 'ciro ' , 

26 US er= ' daniele ' , 

27 password='XXXXXXXXXX' , 

28 workload = 2, 

29 workdir= ' /home/ daniele /requiem/ ' , poll_interval=sleepness /5) 

30 toto=Machine (name= ' toto ' , 

31 maiiager= ' direct ' , 

32 host= ' toto ' , 

33 US er= ' daniele ' , 

34 password='XXXXXXXXXXXX' , 

35 workload = 2, 

36 workdir= ' /home/ daniele /requiem/ ' , poll_interval=sleepness /5) 

37 rosa=Machine (name= ' rosa ' , 

38 manager= ' direct ' , 

39 host= ' rosa ' , 

40 US e r= ' daniele ' , 

41 password='XXXXXXXXXXX' , 

42 workload = 2, 

43 workdir= ' /home/ daniele /requiem/ ' , poll_interval=sleepiiess /5) 

44 mir=Machine (name= ' mir ' , 

45 maiiager= ' direct ' , 

46 host= ' mir ' , 

47 user= ' daniele ' , 

48 password='XXXXXXXXXX' , 

49 workload = 2, 

50 workdir= ' /home/ daniele /requiem/ ' , poll_interval=sleepness /5) 
51 

52 enabled_machines =[f abh , lamh , toto , ciro ,mir , rosa] 
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A. 3 Classi principali: requiem 

A. 3.1 class Job: 

Per ogni punto della griglia viene creato un oggetto Job, la cui funzione principale è quella di 
registrare le proprietà del punto e lo stato del calcolo. Attributi significativi: 

• status -^ Gli stati del calcolo possibili sono: instantiated, ovvero il punto attende di 
essere processato; successf ull, il calcolo è terminato con successo ed i risultati sono stati 
recuperati correttamente; f ailed, il calcolo non è riuscito in quel punto e non verrà più 
tentato. 

• jid -^ Numero unico identificativo 

• name -^ Nome dell'input sul filesystem 

• machine -^ Oggetto Machine a cui è assegnato il lavoro (su quale computer viene svolto) 

• Ipath, Icom, llog, Ichk -^ percorso sul computer locale (dove gira ReQuieM ) in cui si 
trovano o dovrebbero trovarsi, rispettivamente, la cartella contenente la griglia; l'input; il 
log; il checkpoint 

I metodi di questa classe di oggetti sono: 

successo , failureO , notifica il cambiamento di stato e modifica l'attributo status. 

reassign(Machine) , assegna il job ad un oggetto Machine. Nuovi attributi rpath, rlog, rchk, 
rcom definiscono i percorsi dei file sulla macchina remota. 

cmp , questo metodo viene utilizzato da Python per paragonare due oggetti Job e 

stabilire quale è maggiore e quale minore. Gli indici N,M,K sono paragonati. 

A. 3. 2 class JobAgency: 

Questa classe astrae una griglia intera conforme a sez.A: ogni cartella su cui si decide di atti- 
vare la scansione è rappresentata da una JobAgency. Ha metodi di ricerca per restituire un job 
instantiated e attiguo ad uno successfull, e per stabilire quale job venga prima o dopo un altro 
secondo un certo pattern di scansione. La parte più importante della funzionalità di ripristino da 
uno snapshot è codificata in questa classe. 

register(Job) Inserisce l'oggetto Job nella struttura dati (registry) della JobAgency. 

set_workdir("path") Crea una rappresentazione della cartella al percorso path. 

next() Restituisce un oggetto Job eleggibile per il calcolo, oppure una stringa: wait, al momento 
non ci sono punti che soddisfino l'attiguità, ma la griglia non è ancora completa; end, tutti 
i punti della cartella sono completi 

before() Restituisce l'oggetto Job precedente secondo lo schema di scansione. 



A. 3 Classi principali: requiem 

rectifyO Controlla tutti i Job iscritti al registro e riporta a status=instantiated tutti quelli con 
statusT^successful. 

deep() Simile a rectify(). Per ogni Job controlla se esiste il log di Qaussian , e se il calcolo è 
terminato correttamente, e stabilisce in base e questo lo status del Job. 

cmp Gli oggetti JobAgency sono ordinati in base all'attributo priority, ricavato dal 

numero di Job contenuti o imposto manualmente. Le JobAgency con priorità più elevata 
hanno una precedenza parziale nell'occupazione di CPU. 

A. 3. 3 class Machine: 

In questa classe sono presenti tutti i metodi per eseguire lavori su computer remoti, trasferendo 
di volta in volta tutti i file necessari e recuperarando i risultati al termine e nel corso del calcolo. 
Attributi (configurati in ClusterConf ): 

1. manager -^ Indica come devono essere eseguiti i lavori, a scelta tra "direct" e "queue". Se è 
impostato direct, i lavori verranno lanciati in background da shell. L'opzione queue richiede 
invece l'inserimento in una coda gestita dal Grid Engine [48]. 

2. host, port, user, password -^ indirizzo remoto ed autenticazione SSH. 

3. workload -^ numero di CPU/di lavori contemporanei permessi sulla macchina 

4. workdir -^ percorso dove ReQuieM memorizzerà i file temporanei necessari al calcolo 

5. jobs, expired_jobs -^ elenco di lavori in esecuzione / ultimati 
Di seguito i metodi più importanti: 

sshgo(' command") Esegue il comando bash "command" sul computer remoto 

sscp("scp") Trasferisce file tra computer. 

send(Job) Invia il lavoro alla macchina, trasferendo i file con sscp() e lanciando con sshgo() 
comandi conformi all'attributo manager impostato. 

commit(Job) Recupera i risultati dalla macchina e ripulisci il filesystem remoto. Imposta Job.status=failed/s 
in base all'uscita del calcolo. 

monitor(Job) Controlla l'output remoto e stabilisci lo stato di Job. 

poll() Controlla l'output di ciascun lavoro elencato in jobs chiamando monitor (Job). Esegue 
commitMobì ner i lavori ultimati. 



commit(Job) per i lavori ultimati. 
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A. 3.4 class Requiem: 

Questa classe implementa i punti 111,1-9 dello schema di esecuzione presentato all'inizio dell'ap- 
pendice. 

update() Ricorstuisce le statistiche e le stampa sul file panel.log: 

### [ReQuieM] INFORMATION PANEL [ReQuieM] ### 
Wed Mar 19 13:11:29 2008 



CONFIGURED MACHINES: 



NAME MANAGER 
f abh queue 
lamh queue 
toto direct 
Cluster Capacity: 12 



HOST USER 

hydra fabrizio 

hydra lami 

toto daniele 



WL CONNECTION 
5 OK 



JOB AGENCIES: 

NAME 
0002 
0001 



WORKDIR PRIORITY MODE JOBS COMPLETION 

/scr/daniele/topist6/000 10 Isp {1:19,2:95,3:1995} 94.59% 
/scr/daniele/topist6/010 9 Isp {1:19,2:95,3:1995} 94.55% 



SUBMITTED JOBS: 



ERRONEOUS JOBS: 



STATUS MACHINE AGENCY 



STATUS MACHINE AGENCY LAST LINE 



refresh_agencies() Apre un file di comunicazione e crea/rimuove JobAgency rispondendo ai 
comandi letti. 

place() Ricerca una Machine libera un nuovo lavoro. 

cycle() Questo metodo è il metodo chiamato iterativamente da RequiemDirector . Aggiorna 
le JobAgency con refresh_agencies() (riga:2). Trova lavori da svolgere interrogando tutte 
le JobAgency, riga 6 - l'iteratore for richiama ad ogni ciclo il metodo next() della classe 
JobAgency. Se la risposta della JobAgency è un oggetto Job, esso viene assegnato ad 
una maccina dal metodo place(); se è una stringa tra wait oppure end, viene esaminata la 
JobAgency successiva. Al termine dell'esecuzione cycle() aggiorna le statistiche di scansione 
con update(). Nel caso nessuna macchina abbia slot disponibili, richiede l'aggiornamento 
dei tutte le macchine e dello stato dei lavori in esecuzione. 

1 def cycle(self ): 



2 


self . refresh_agencies () 


3 


if len(self . agencies) = = 0: 


4 


return False 


5 


waiting=0 


6 


for AG in self . agencies : 


7 


placed=True 


8 


for job in AG : 


9 


if repr ( type ( j ob)) = ="<type 'instance'>" 


10 


placed=False 


11 


placed=self .place(job) 



90 



A. 4 Esecuzione da riga di comando: RequiemDirector 

12 if placed=False : 

13 AG . log ( ' Waiting a free machine ', 1 = 1) 

14 break 

15 [ 1 

16 elif job='wait': 

17 AG . log (' Waiting for parent jobs to end ',1 = 1) 

18 waitiiig+=l 

19 break 

20 elif job=='end': 

21 AG.destroyO 

22 self . agencies . remove ( AG) 

23 self . empty_agencies+=[AG ] 

24 break 

25 if placed^=False : 

26 break 
27 

28 if waiting=len ( self . agencies ) : 

29 self . update (mac=True ) 

30 else : 

31 self . update () 
32 

33 return True 

rectify(), deep_resume() Esegue le funzione omonime già viste nella classe JobAgency, ma su 
tutte le JobAgency configurate. 

A. 3. 5 class LogFrame: 

Utilità per il logging interno di ReQuieM , che produce un registro con vari livelli di guardia in 
main_logfile, del tipo: 



[0]@Tue Mar 18 15:30:20 2008 
[2]@Tue Mar 18 15:31:40 2008 
[l]@Tue Mar 18 15:31:40 2008 



SSH connection was lost [m:fabh] 

Job enqueued at fabh [j:j0652] 

Job Object Reassigned to fabh [j:j0671] 
VVVV agency=0001 jid=j0671 

IMI lpath=/scr/daniele/topist6/010/ Ilog=/scr/daniele/topist6/010/x5-y4-zl2.1og 

IMI Icoin=/scr/daniele/topist6/010/x5-y4-zl2.com Ichk=/scr/daniele/topist6/010/x5-y4-zl2. chk .bz2 

IMI rpath=/home/f abrizio/daniele/requiein/j0671/ rcoin=/home/f abrizio/daniele/requiem/j0671/x5-y4-zl2 .com 

>>>> rlog=/hoine/f abrizio/daniele/requiein/j0671/x5-y4-zl2.1og rc]ik:=/hoine/f abrizio/daniele/requiem/j0671/x5-y4-zl2.chk .bz2 [j : J0671] 



[2](§Tue Mar 18 15:32:22 2008 

[l]@Tue Mar 18 15:32:22 2008 

[2]@Tue Mar 18 15:32:42 2008 

[2]{§Tue Mar 18 15:33:03 2008 

[l](§Tue Mar 18 15:33:03 2008 

[l](§Tue Mar 18 15:33:03 2008 

[l]@Tue Mar 18 15:33:03 2008 



Job enqueued at fabh [j:j0671] 
Waiting for parent jobs to end [ag:0001] 
Job started by remote queue [j:j0652] 
Job started by remote queue [j:j0671] 
Poli: 2 running, committed jobs [m:fabh] 
Poli: running, committed jobs [m:lamh] 
Poli: running, committed jobs [m:toto] 



A. 4 Esecuzione da riga di comando: RequiemDirector 

RequiemDirector ^ è il programma eseguibile che può essere lanciato da riga di comando per 
avviare la scansione. Senza opzioni a riga di comando ReQuieM ricostruirà da zero tutte le sue 
strutture dati, ricominciando daccapo la scansione. 

L'opzione "resumé" ripristina la scansione dal punto in cui il programma è uscito l'ultima volta. 
Le opzioni "rectify" e "resumé" sono utili per correggere le inconsitenze tra registro dei lavori e 
calcoli effetivamente svolti che possono emergere a seguito di una interruzione impovvisa. 



"^Director, direttore d'orchestra, con riferimento all'origine musicale del nome ReQuieM , poiché il compito di 
questo file è governare l'avvio, l'esecuzione e la chiusura del programma. 
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L'interruzione del programma può infatti lasciare dei Job in uno stato indeterminato, inducendo 
ReQuieM ad ignorarli. Oppure può fermare a metà il trasferimento di file. Con queste due 
opzioni si può riparare lo stato degli oggetti. 

Una volta avviato ReQuieM in base a queste opzioni, RequiemDirector esegue sostanzial- 
mente questo ciclo: 

1 cy=True 

2 while cy : 

3 print 'CYCLE' 

4 save ( RequiemObj ) 

5 t l=t ime ( ) 

6 try : 

7 cy=RequiemObj . cycle 

8 t2=time() 

9 w=l+sleepness /(t2-tl) 

10 if w>sleepness : 

11 w=sleepness 

12 sleep (w) 

13 except : 

14 cy=False 

15 e=f ormat_exc 

16 RequiemObj . log( 'EXCEPTION RAISED IN MAIN MANAGER: \ n %s ' % e ,1=3) 

17 save ( RequiemObj , main_dump+ ' . emergency ' ) 

18 print 'Emergency core dump finished... ' 

Alla riga 3 viene salvata la gerarchia di oggetti dell'istanza della classe Requiem per la sessione 
corrente, ovvero l'oggetto RequiemObj, nel file main_dump. 

In riga 7 viene chiamato il metodo RequiemObj. cycle() che, se possibile, propaga di uno step 
la scansione. Quindi RequiemDirector attende il tempo necessario a far passare l'intervallo di 
tempo sleepness (A. 2) dall'inizio del ciclo. 

Uno dei punti di forza di Python è l'estrema flessibilità e tolleranza nel gestire le eccezioni. Il 
costrutto try: ... / except: ... permette salvare un ultimo snapshot immediatamente dopo che 
l'eccezione è stata sollevata. Il paragone tra main_dump e main_dump. emergency può quindi 
aiutare a determinare meglio le circostanze di crash. 

Il programma ReQuieM è stato di fondamentale importanza per acquisire in tempi utili i dati 
presentati in questa tesi. 
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Il programma GeoM<n è stato sviluppato per manipolare la OJtairice Z di un grande numero di 
geometrie molecolari e per movimenti che accoppiano in modo complesso molteplici coordinate 
interne. Il nome GeoM<n è infatti una abbreviazione per Geometry Manipulator. 

Di seguito è descritta brevemente l'implementazione nel linguaggio di programmazione Python 
[39], e l'organizzazione interna in funzioni, classi di oggetti e metodi. Quindi sono presentati 
alcuni semplici script che utilizzano GeoM<n per risolvere le difficoltà incontrate nel corso del 
lavoro di tesi. 

B.l class Mol: 

L'implementazione consiste in una classe di oggetti Mol ed un numero di funzioni accessorie. 
La classe Mol definisce attributi per memorizzare la OJlaitice Z e le coordinate cartesiane delle 
molecole; i metodi della classe permettono poi alcune trasformazioni. 

Attributi principali 

• types -^ Lista contenete i tipi di atomi, nell'ordine delle coordinate cartesiane 

• coords -^ Lista con le coordinate cartesiane degli atomi 

• zm -^ Lista con le righe della OJtairice Z simbolica, [[elementi riga l],[elem. riga 2],etc] 

• var -^ Dizionario delle variabili della OJlaitice Z simbolica, nomeivalore. 

• zmn -^ OJlairice Z numerica 

Metodi 

loadc("path"), loadz("path") Legge un file contenente la struttura molecolare dal percorso "pa- 
th" trasferendone l'informazione in strutture dati Python. loadc() interpreta il file come 
coordinate cartesiane; loadz() come 9}tairice Z . 

internal(A,B7C,D) Dati A,B^(7;/? vettori di punti in coordinate cartesiane, il metodo restituisce 
le coordinate interne tra di loro. Almeno due punti vanno dati come input, ed in questo caso 
la funzione restituisce solo la distanza tra di loro. Per tre punti ritorna distanza e angolo; 
per quattro, anche il diedro. 
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seq2zm(seq) Genera una OJtatrice Z dalle coordinate cartesiane precedentemente caricate con 
loadc(). Ogni atomo della OJlatrice Z è definito rispetto ai tre che lo precedono. Il parametro 
opzionale seq è una lista che riporta l'ordine con cui devono essere disposti gli atomi nella 
matrice. Ad esempio, se le coordinate cartesiane sono 4, e si desidera una 9}tairice Z con il 
primo e l'ultimo atomo scambiati, seq=[4,2,3,l]. 

fill(zm) Date le coordinate cartesiane degli atomi ed un modello di completo di 2Jlatrice Z per 
la molecola, calcola il valore di ogni variabile in OJlatrice Z in modo da corrispondere alla 
geometria espressa dalle coordinate cartesiane. Questo permette di generare 9Jlatrici Z 
complesse, in cui le definizioni di atomi non sono date esclusivamente tra atomi successivi. 

corder(dict) Riordina le coordinate cartesiane secondo il dizionario dict: {nuova posizione: vec- 
chia posizione}. Ad esempio, volendo scambiare il primo e l'ultimo di 4 atomi, dict={l:4,2:2,3:3,4:l}. 

zm2cart( format) Genera una rappresentazione in coordinate cartesiane da quella in OJtaitice Z 
. Il parametro format stabilisce se la OJtattice Z è in formato numerico o letterale. 

L'inizializzazione di un oggetto Mol può richiedere direttamente il caricamento di una geometria: 
obj=Mol(zmat="geo.com",xyz="geo.xyz"). 

B.2 Funzioni accessorie 

leng(a), normalize(a), angle(a,b,c), dihedral(a,b,c,d) 

La prima calcola la lunghezza del vettore a, la seconda lo normalizza. angle() e dihedral() Ri- 
chiedono come argomenti le posizioni cartesiane degli atomi e restituiscono la coordinata interna 
indicata dal nome della funzione. 

locate(A,B,C,r,a,d), arbrot(P,A,a) 

locate() restituisce un vettore con distanza r da A; angolo a con A,B; diedro d con A,B,C. 

arbrot(P,A,a) ruota dell'angolo a il punto P attorno all'asse passante per l'origine e Ax ed è 
utilizzata da locate per calcolare il diedro, utilizzando il metodo dei quaterioni [14] 

zrender(Mol), cartrender(coords, types) 

Stampano la OJlaince Z (zrender) e le coordinate cartesiane (cartrender) in formati direttamente 
utilizzabili da Qaussian o da programmi di visualizzazione, come Molden. 

sa ng3 ( A , B , C, D , ri fi , rif2, die, r) 

Questa funzione implementa il metodo per il posizionamento degli atomi fittizi esposto in sez.2.1.2. 

1 def sang3 (A , B , C ,D , rif l=None , rif 2=None , die=None , r = 1.) : 

2 """Determine the position and angle of a point X that has 

3 distance r from D and the same angle X-D-(A,B,C) . " " " 
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B.3 Caso d^uso: aggiunta degli atomi fittizi 

4 A=iiormalize (D— A ) ; B=iiormalize (D— B ) ; C=iiormal ize (D— C ) ; 

5 S=array ( [ A , 

6 B, 

7 C ]) 

8 K=array([ [1.] , 

9 [1.], 

10 [1.] ]) 

11 try: 

12 iS=inv(S) 

13 X=dot(iS,K) 

14 X=array([ X[0][0], X[1][0], X[2][0] ]) 

15 a=arccos ( 1 . / leng (X ) ) 

16 if die ! = None : 

17 r i f l=D-r i f 1 ; r i f 2=D-r i f 2 

18 besta=360 

19 rd=dihedral(X,0 ,rif 1 ,rif2) 

20 if rd*die<0: 

21 X=-X 

22 X=X*r*cos(a) 

23 except LinAlgError : 

24 #Three vectors are on the same piane . . . very simple ! 

25 print 'LinAlgError' 

26 a=pi/2 

27 X=normalize ( cross ( A , B ) ) * r 

28 X=X+D 

29 return X , a 

Dati l'atomo interanello D e gli atomi vicinali A,B,C, come coordinate cartesiane, sangS calcola 
la posizione dell'atomo X avente distanza r da D e angoli XD{A^B^C) uguali tra di loro. Gli 
argomenti opzionali rifl, rif2 sono atomi rispetto ai quali X deve formare il diedro die XDriflrif2 

B.3 Caso d'uso: aggiunta degli atomi fittizi 

La prima applicazione di GeoM<n è stata la creazione della OJlatrice Z progettata in sez.2.1.4 
per tutti i punti a disposizione in coordinate cartesiane (tra cui quelli elencati in tab.2.3). Per lo 
scopo abbiamo scritto uno script: 

1 from geoman import * 

2 def dummif y ( xyz= ' minEsO . xyz ' , zmat= ' minEsO . com ' ) : 

3 """Add dummy atoms as needed""" 
4 

5 m=Mol ( xyz=xyz , zmat=zmat ) ^ Load both cartesian and internai coordinates 

6 T^Calculate dummy atoms positions and add them to cartesian coords 

7 # The first dummy must bave the same angle with 3,5,6 and dist 1 from 2: 

8 X , a=sang3 (m. coords [3] ,m. coords [16] ,m. coords [13] ,m. coords [12]) 

9 m . coords+=[X ] 
10 m. types+=[ '0 ' ] 
11 

12 #The second dummy must bave the same angle with 2,7,8; 

13 #dist 1 from 3; 

14 #plus have correct dihedral angle from 2 and the former dummy: 

15 X , a=sang3 ( m. coords [12] ,m. coords [4] ,m. coords [2] ,m. coords [3] ,\ 

16 m.coords[12],m.coords[23],m.var[ 'DX' ] ) 

17 m. coords+=[X] 

18 m. types+=[ '0 ' ] 
19 
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B GeoM<n 

20 

21 m . corder ( z2xyzMolCas ) 7^ reorder cartesian coordinates in the same order zmatrix has 

22 7^. corder ( z2xyz ) 

23 Izm , Ivar , nzm= m.fill() ^ Fili the template with info from cartesian coords 

24 m.zm=lzm; m.var=lvar; m . zmn=nzm ^ And memorize fi lied zmatrix 

25 zrender (mol=m , f ormat= ' molden ' ) 
26 

27 f or geom in (xyzl , xyz2 ,...): 

28 dummif y ( xyz=geom [0] , zmat=" de fault . zmat . com" ) 

In riga 2 si definisce la funzione che sarà poi applicata in riga 24 a tutti i punti per cui si desidera 
creare la OJtaitice Z , dummify. Questa inizializza l'oggetto Mol m (r.5), aggiunge gli atomi dummy 
con la funzione sangS (r.8,15) alle coordinate cartesiane. Infine viene compilata la filatrice Z sulla 
base dello schema default.zmat.com definito in sez.2.1.4 (r. 18-21) e stampa il risultato (r.22) 

B.4 Caso d'uso: calcolo della distanza tra gli H dei due anelli 

La figura 2.10 è riporta le distanze tra gli atomi di idrogeno tra i due anelli, calcolate per molte 
geometrie molecolari con lo script seguente in modo molto più rapido rispetto alla misurazione 
"manuale" con programmi di visualizzazione come Molden o VMD. 

1 from geoman import * 

2 m=Mol ( zmat=" ci — diepir . com" ) 

3 hea='X 14-19 14-23 14-24 13-19 13-23 13-24 18-19 18-23 18-24' 

4 print hea 

5 ent='%i\t '+'\t ' . join([ '%5.3f ']*9) 
6 

7 for da in range ( , 360 , 1 1 ) : 

8 m. var [ 'D3'] = d3 

9 m . zm2cart () 

10 print ent % (da ,m. internai (14 ,19) [0] , m . internai ( 1 4 , 2 3 ) [0] , \ 

11 m. internai (14,24)[0], m. internai (13,19)[0], m. internai (13,23)[0], \ 

12 m. internai (13,24)[0], m. internai (18,19)[0], m. internai (18,23)[0],\ 

13 m. internai (18 ,24) [0]) 
14 

15 m=Mol ( zmat=" ci — diepir 2 .com") 

16 print '\n\n\n'+hea 

17 for da in range ( , 360 , 1 1 ) : 

18 m. var [ 'D3'] = da 

19 m . zm2cart () 

20 print ent %(da , m . internai ( 1 4 , 1 9 ) [0] , m . int ernal ( 1 4 , 2 3 ) [0] , \ 

21 m. internai (14,24)[0], m. internai (13,19)[0],m. internai (13,23)[0], \ 

22 m. internai (13,24)[0], m. internai (18,19)[0], m. internai (18,23)[0],\ 

23 m. internai (18 ,24) [0]) 

I file "ci-diepir.com" e "ci-diepir2.com" contengono 2 OJlaitici Z per geometrie con P3 distorto 

in modo diverso. Queste vengono caricate nelle righe 2,12. 

II programma produce questo tipo di output per ciascuna geometria di partenza: 



X 


14-19 


14-23 


14-24 


13-19 


13-23 


13-24 


18-19 


18-23 


18-24 





4.429 


1.768 


2.592 


5.037 


3.307 


3.342 


2.259 


4.689 


4.828 


11 


4.506 


1.780 


2.336 


5.028 


3.410 


3.260 


2.362 


4.568 


4.859 


22 


4.548 


1.878 


2.106 


4.990 


3.543 


3.217 


2.525 


4.418 


4.855 


33 


4.553 


2.046 


1.919 


4.925 


3.699 


3.218 


2.732 


4.242 


4.817 


44 


4.522 


2.265 


1.796 


4.833 


3.868 


3.262 


2.968 


4.042 


4.746 


55 


4.455 


2.514 


1.758 


4.716 


4.044 


3.346 


3.217 


3.824 


4.642 


352 


4.350 


1.816 


2.785 


5.025 


3.255 


3.424 


2.229 


4.757 


4.785 
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Elenco delle figure 



1.1 L'occhio umano 2 

1.2 A sinistra: Stratificazione della retina (luce dal basso). A destra: Disposizione controintuitiva degli strati 
retinici, con le cellule sensibili a contatto con l'epitelio retinico pigmentato 3 

1.3 Schema di una cellula bastoncello [5]. La luce proviene dall'alto 4 

1.4 Rappresentazione della Rodopsina transmembrana, in rosso il cromoforo. A destra: legame imminico (base 

di Schiff protonata) con la lisina. [8] 5 

1.5 Isomerizzazione del retinale nella cavità opsinica, incapace di legare la forma ali-trans [40]. Conformazioni 
intermedie della Rodopsina prima di espellere completamente il retinale [21]. La Fotorodopsina è una 
forma calda di Batorodopsina, in cui il retinale è altamente distorto. La sua esistenza è tuttora discussa . 6 

1.6 Cascata biochimica che sfocia nell'iperpolarizzazione della membrana cellulare del bastoncello. Step): 
Fotoreazione da Rodopsina a Metarodopsina. 2) Attivazione della Proteina G. 3) Produzione del complesso 
inibitorio tra una subunutà a di G ed il GTP. 4) Inibizione della fosfodiesterasi PDF. 5) Il cGMP viene 
esaurito dalla PDF ed i canali sodio non possono più rimanere aperti. [28] 6 

1.7 Spettroscopia laser al femtosecondo con metodo pump-probe (pompa-sonda) su un getto molecolare. [64] 7 

1.8 A sinistra: stato eccitato e fondamentale per la molecola NaI. A destra: Oscillazioni coerenti della distanza 
internucleare Na-I sulla superficie di potenziale dello stato eccitato. [64] 8 

1.9 Sinistra: cambiamento di assorbimento (6T/T) a 4 lunghezze d'onda. Destra: spettri differenziali in 4 
tempi. [47] 9 

1.10 A sinistra: Spettro Raman a tempi diversi. A destra: Modello dinamico di Kukura [25] 10 

1.11 Cammini di energia minima sullo stato eccitato per il retinale [13] 11 

1.12 A sinistra: strutture complesse in grado di compiere semplici movimenti quando viene attivato uno switch 
presente sulla catena o sugli anelli [1]. A destra: braccio molecolare con estensione variabile basato su un 
diarilidene ed implementazione in un controllo di giunzione a livello nanometrico [9] 12 

1.13 Diagramma del lavoro di Tesi 15 

2.1 La nostra numerazione del CPP (idrogeno e dummies esclusi). Conformazione trans. Gli atomi 3,2 sono 

detti interanello o giunzione 16 

2.2 Sinistra: asse torsionale del Z-penta-3,5-dieniminio [36] / Destra: Coordinata di piramidalizzazione per il 
4-Me(C6H9)CHF [24] 17 

2.3 Nomenclatura dei legami sul sistema tt 18 

2.4 Definizione dell'atomo fittizio e della coordinata P 21 

2.5 Proiezioni di Newman raffiguranti l'angolo di torsione. Da sinistra: geometria trans, CI, angolo X . . . . 24 

2.6 I primi 8 atomi definiti nella OJtatrice Z contengono la coordinata di torsione e le due di piramidalizzazione 25 

2.7 Distorsione sul legame 9-10 provocata da P2 27 
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2.8 A sinistra: rappresentazione geometrica delle tre variabili che determinano la posizione degli atomi 5 e 6. 

A destra: sistema cartesiano per il ricalcolo degli angoli diedri esterni 28 

2.9 Geometria ottenuta ponendo P2 = 110° e X=209° e applicando la correzione degli angoli esterni. Gli atomi 
d'idrogeno 14-19 sono molto ravvicinati 30 
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